Investigation on Thermal-Hydraulic Performance in a Printed Circuit Heat Exchanger with Airfoil and Vortex Generator Fins for Supercritical Liquefied Natural Gas

Abstract

In this study, the thermal and hydraulic performances of six different printed circuit heat exchanger (PCHE) configurations with supercritical liquefied natural gas as the working fluid are studied by a numerical method. Firstly, the thermal and hydraulic performances of the PCHE with airfoil fins are investigated at different operating pressures, which indicate that the PCHE with airfoil fins has better thermal performance but worse hydraulic performance when it operates at higher pressure condition. Furthermore, the effects of different PCHE configurations on the thermal and hydraulic performances are analyzed in detail. The results show that delta vortex generators (VGs) arranged upstream the airfoils can provide better thermal performance than VGs arranged downstream the airfoils, and the thermal performance of the common-flow-up configuration is not always better than the common-flow-down configuration. Finally, the overall heat transfer performance of PCHEs is evaluated by the volume goodness criterion under the same mass flow rate. The results show that the straight channel PCHE with front common-flow-up delta VG and airfoil fins has the best volume goodness factor among the six PCHE configurations.

Introduction

Natural gas (NG) is well-known as one of the cleanest form of fossil fuels, and it is the fastest growing primary energy source owing to its advantages such as clean and efficient combustion, reliable and durable supply, energy and specific density, convenient and clean usage, and flexible implementations (vehicles, power generation plants, industrial and residential usage, etc.) [1]. The liquefied natural gas (LNG) technology provides an economically feasible way of transporting NG over long distances and currently accounts for nearly 30% of the international trade of this resource [2]. Floating liquefied natural gas is an alternative LNG supply method, and it is easy to produce, carry and import of the LNG as it is compared to the traditional LNG supply chain for the offshore reserves [3], and liquefied natural gas - floating production storage and offloading (LNG-FPSO) is a promising energy-efficient method compared to onshore NG treatment facility when the gas source is located in a distant offshore ocean [4]. A further development is the floating regasification units that transform LNG back to NG at the market location, and such units are called liquefied natural gas - floating storage and regasification units (LNG-FSRU). Process technology for LNG-FSRU can be presented as follows: LNG is transferred into cargo tank by dedicated LNG feed pump, then LNG is sent to high pressure booster pump which is used to make high pressure before entering LNG vaporizer, after that LNG will be vaporized and sent to end user [5]. Four main kinds of LNG vaporizers are commercially available, that is, intermediate fluid vaporizer, open rack vaporizer, super open rack vaporizer, and submerged combustion vaporizer. However, the heat exchangers for LNG-FPSO and LNG-FSRU should be more compact than the ones in onshore liquefaction plant due to its limited platform area [4]. Printed circuit heat exchanger (PCHE) is a newly developed heat exchanger with favorable characteristic of compactness, and its potential to further penetrate the market is rather substantial, especially if the heat exchanger size, weight and endurance are major concerns [2]. Therefore, PCHE could be used for LNG-FRSU as a LNG vaporizer.

The first PCHE was originally developed in 1980 at the University of Sydney, while a major PCHE manufacturer and pioneer is the Heatric company (UK), which is highly specialized in the area since 1985 [6]. PCHEs are manufactured by photochemical etching [7] and diffusion bonding technologies, and are able to reliably operate under large temperature range (from −200 °C to 900 °C) and high pressure (60 MPa) conditions [8]. During the past few years, researchers have spent considerable effort to study the thermal and hydraulic performance of PCHEs. The common flow channels of PCHEs are straight channels, wavy channels, zigzag channels, S-shaped channels and channels with airfoils. Table 1 lists references about the thermal and hydraulic performance of PCHEs with different configurations [9–23]. From the above review, it is concluded that compared with the straight channels, the heat transfer performance may be enhanced by non-straight channels, especially channels with airfoil fins. But most of references in Table 1 used supercritical CO₂ and helium as working fluids, and there are few studies applying PCHE in LNG regasification process, it is necessary to study the thermal and hydraulic performance of airfoil fin PCHEs with supercritical LNG as the working fluid in supercritical LNG heating process.

Table 1. List of references about different PCHE channels.

References	Channels	Working fluid	Working conditions	Main conclusions
Mylavarapu et al. [9]	straight	helium	cold side: 358–663 K and 1–2.7 MPa; hot side: 481–1063 K and 1–2.7 MPa; mass flow rate: 15–49 kg/h	The critical Reynolds number from the laminar to transition flow regime was about 1700, and the experimental correlations were given.
Chu et al. [10]	straight	supercritical CO2, water	hot supercritical CO2: 310–375 K, 8–11 MPa, and 150–650 kg/h	PCHE had better comprehensive performance at higher pressure, and the comprehensive performance was significantly reduced nearly 17.6% at the transcritical state.
Khan et al. [11]	wavy	helium	cold side: 813 K and 3 MPa; hot side: 1173 K and 3 MPa; Re: 350–2100	The wavy channel could provide better thermal-hydraulic performance compared with the straight channel.
Aneesh et al. [12]	wavy	helium	cold side: 813 K, 3 MPa, and Re: 451–2256; hot side: 1173 K, 3 MPa, and Re:348–1741	The highest values of thermal advantages of trapezoidal, sinusoidal and triangular wavy channel, compared with the straight channel, were 41, 33, and 28%, respectively.
Lee and Kim [13]	zigzag	supercritical CO2	cold side: 396 K, 8.3 MPa, and Re: 67000–280000; hot side: 411.2 K, 2.5 MPa, and Re: 71500	The effectiveness of the heat exchanger was maximized when the cold channel angle was similar to the hot channel angle.
Ma et al. [14]	zigzag	helium	cold side: 813 K, 3 MPa; hot side: 1173 K, 3 MPa; mass flow rate: 0.072–0.324 kg/h	With the increase of inclination angle, the heat transfer and pressure drop both increased.
Pan et al. [15]	zigzag	supercritical LNG, propane	cold side: 128.9 K, 7.3 MPa, and 0.0002995 kg/s; hot side: 293.15 K, 0.64 MPa, and 0.00045 kg/s	The larger the bending angle was, the better the heat transfer performance was accompanying larger pressure drop. The channel bending angle 15° offered the best comprehensive performance.
Ngo et al. [16]	S-shape	CO2 , water	water: 280 K, 0.25 MPa, and 48 kg/h; CO2: 391 K, 11.5 MPa, and 57.8 kg/h	The S-shaped PCHE provided about 3.3 times less volume, and lower pressure drop by 37% in the CO2 side and by 10 times in the H2O side, compared with the existing supplier.
Tsuzuki et al. [17]	S-shape, zigzag	supercritical CO2	cold side: 363–381 K and 6.5–10 MPa; hot side: 553–573 K and 2.2–3.2 MPa; mass flow rate: 40–80 kg/h	The S-shaped PCHE, compared with a conventional zigzag configuration, had one-fifth of the pressure drop with equal thermal-hydraulic performance.
Xu et al. [18]	airfoil	supercritical CO2	inlet temperature: 303 K, inlet pressure: 8 MPa, and wall temperature: 333 K	With increase in mass flow rate, the overall thermal and hydraulic performance of airfoil fin PCHE was better than that of other types of fin.
Ma et al. [19]	airfoil	supercritical CO2	cold side: 761.96 K and 7.83 MPa; hot side: 1173.15 K and 7.59 MPa	The maximum difference of Nu between the two models with and without airfoil fin in the PCHE was up to 6.7%, and the maximum difference of f-factor was 6.4%.
Chu et al. [20]	airfoil	supercritical CO2	inlet pressure: 8 MPa, and wall temperature: 393 K	The heat transfer rate decreased because of strong variation of properties of supercritical CO2, but the pressure loss remained almost constant along the main flow direction.
Cui et al. [21]	airfoil	supercritical CO2	cold side: 381.05 K, 8.28 MPa, and 0.4728–1.8912 g/s; hot side: 553.05 K, 2.52 MPa, and 0.4335 g/s	The staggered arrangement and the proper shape of airfoil fin could reduce the effect of boundary layer and improve the heat transfer.
Chu et al. [22]	airfoil	supercritical CO2	inlet temperature: 473.15 K, inlet pressure: 8 MPa, and wall temperature: 313 K	Nu of airfoil fins with different cambers could be improved by 42, 54, 75, and 66%, besides f-factor averagely increased 125, 460, 980. and 890%, respectively.
Zhao et al. [23]	airfoil	supercritical LNG	inlet temperature: 121 K, inlet pressure: 10.5 MPa, and heat flux: 58713.75 W/m^2	A sparser staggered arrangement of fins showed a better thermal and hydraulic performance in airfoil fin PCHE.

On the other hand, vortex generators (VGs) can be mounted on the channel wall to continuously disturb the viscous sublayer along the main flow direction and to produce longitudinal vortices (LVs) inducing strong swirling motion, therefore, it may greatly enhance heat transfer [24]. Following Mehta et al. [25], LVs generated by VGs can be denoted as common-flow-down configuration and common-flow-up configuration, respectively. Wu et al. [26] numerically investigated the thermal and hydraulic performance of a heat exchanger with common-flow-up interrupted delta winglet VGs. The results showed that compared with a heat exchanger without VGs, the heat transfer capacity and overall performance increased by 35–60% and 19–64%, respectively. Naik and Tiwari [27] numerically investigated effects of various geometrical parameters of common-flow-down rectangular winglet VGs on heat transfer enhancement in flow over a flat plate. The results indicated that increase in enhancement in heat transfer was observed with increase in length of VGs, and the overall heat transfer performance was the best at a particular attack angle of VGs. A comparative study of effects of elliptical pole parameters, the thickness, the length, the height, the attack angle, the transverse pitch and the longitudinal pitch of common-flow-up winglet VGs on the thermal and hydraulic performance in a rectangular channel was numerically investigated by Tang et al. [28]. The results showed that the intensity of heat transfer could be greatly increased by the increase of length, height, attack angle and transverse pitch of VGs, accompanied with an increase of pressure drop. Tang et al. [29, 30] compared the effect on the thermal and hydraulic performance of common-flow-down winglet VGs and common-flow-up winglet VGs. The results showed that the thermal performance and hydraulic performance of common-flow-up VGs was better than that of common-flow-down VGs. From the above review, it was found that these two configurations both greatly improved the heat transfer enhancement. Therefore, the combination of VGs and airfoils may bring greater heat transfer enhancement. In this study, six PCHE models are used to study the thermal and hydraulic performance in the supercritical LNG heat exchange process. The results may provide a practical guidance on heat transfer enhancement method for LNG vaporization devices.

Numerical method

Physical model

Methane is the major chemical component of LNG which varies between 87.0 and 99.8% [31] for different LNG sources. The composition of LNG in China is shown in Table 2 [32]. The mean relative deviations of thermal-physical properties of LNG and methane, that is, specific heat c_p, density ρ, thermal conductivity λ and dynamic viscosity μ, in the studied temperature range at 8 MPa are 9.0, 8.6, −4.8, and 5.7%, respectively, and it is acceptable for engineering computation and analysis. Therefore, LNG can be assumed as 100% methane [33–35], and the thermal-physical properties of methane are used for LNG. As is well-known, gas exists as a single phase medium when the pressure and temperature are above the critical state, so the phase change phenomenon of methane no longer happens when its pressure and temperature is above 4.59 MPa and 190.56 K. The thermal properties of methane, that is, specific heat c_p, density ρ, thermal conductivity λ and dynamic viscosity μ, are calculated by REFPROP at different pressures. In comparison with the traditional fluid flow and heat transfer behavior, the thermal-physical properties of supercritical methane change dramatically with temperature and pressure as shown in Figure 1. The thermal-physical properties variations with temperature can be expressed by the user defined function (UDF) in Ansys Fluent 19.1, and the errors between data from REFPROP and the computed data by UDF are within ±2%.

Figure 1. Thermal-physical properties of methane at different pressures.

Table 2. Composition of LNG [32].

Composition	Methane	Nitrogen	Ethylene	Propane	Butane
Fraction %	89.62	3.05	5.83	0.14	1.36

A three-dimensional heat exchanger unit is designed to implement the numerical calculations, as shown in Figure 2. The methane-flow is described in a coordinate system, in which x, y, z are streamwise, spanwise, and normal coordinates, respectively. The channel height H is 1.0 mm, and the channel width W is 2.0 mm. In order to support a uniform velocity at the inlet and suppress the backflow at the outlet, the inlet and the outlet are both extended to 16 times of the channel height. Therefore, the channel lengths (L₁, L₂, L₃) are 16.0 mm, 64.0 mm and 16.0 mm, respectively.

Figure 2. Schematic diagram of the physical model.

The NACA 23024 airfoil fin is used in the present study, as shown in Figure 3. The geometrical parameters of the airfoil fin, that is, the maximum thickness W_f, the chord length L_f and the height H_f, are 0.96, 4.0, and 1.0 mm, respectively.

Figure 3. Schematic diagram of the airfoil.

Six different PCHEs are compared: Case A, straight channel PCHE; Case B, straight channel PCHE with airfoil fins; Case C, straight channel PCHE with front common-flow-up delta VG and airfoil fins; Case D, straight channel PCHE with front common-flow-down delta VG and airfoil fins; Case E, straight channel PCHE with rear common-flow-up delta VG and airfoil fins; Case F, straight channel PCHE with rear common-flow-down delta VG and airfoil fins. The list of different configurations is shown in Table 3, and the geometric parameters are listed in Table 4.

Table 3. List of different PCHE configurations.

PCHE	Characteristics	Schematic Diagram
Case A	straight channel
Case B	airfoil
Case C	front common-flow-up delta VG and airfoil
Case D	front common-flow-down delta VG and airfoil
Case E	rear common-flow-up delta VG and airfoil
Case F	rear common-flow-down delta VG and airfoil

Table 4. Geometric parameters of different PCHE configurations.

Parameter	Symbol/unit	Case A	Case B	Case C	Case D	Case E	Case F
Channel length	L2/mm	64.0	64.0	64.0	64.0	64.0	64.0
Channel width	W/mm	2.0	2.0	2.0	2.0	2.0	2.0
Channel heigth	H/mm	1.0	1.0	1.0	1.0	1.0	1.0
Transverse pitch	Pt/mm	–	–	1.0	1.0	1.0	1.0
Longitudinal pitch	Pl/mm	–	8.0	8.0	8.0	8.0	8.0
Max airfoil thickness	Wf/mm	–	0.96	0.96	0.96	0.96	0.96
Airfoil chord length	Lf/mm	–	4.0	4.0	4.0	4.0	4.0
Airfoil heigth	Hf/mm	–	1.0	1.0	1.0	1.0	1.0
VG attack angle	α/deg	–	–	45	135	45	135
VG length	LVG/mm	–	–	0.5	0.5	0.5	0.5
VG thickness	WVG/mm	–	–	0.05	0.05	0.05	0.05
VG height	HVG/mm	–	–	0.8	0.8	0.8	0.8
Row	N	–	8	8	8	8	8

The governing equations for continuity, momentum and energy may be expressed in tensor notation as (1) $$\nabla \cdot \left(\rho \overrightarrow{U}\varphi \right)=\nabla \cdot \left({\Gamma }_{\varphi }\nabla \varphi \right)+S_{\varphi }$$(1)

In the above equation, the dependent variable, ϕ, stands for the velocity components, temperature, k and ε. The terms Γ_ϕ and S_ϕ represent the appropriate diffusion coefficients and the source terms, respectively. The particular expressions for ϕ, Γ_ϕ and S_ϕ [36] are summarized in Table 5.

Table 5. Expressions for ϕ, Γ_ϕ, and S_ϕ.

ϕ	Γϕ	Sϕ
1	0	0
u	μ + μt	$S_{\mathrm{u}}=-\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}((\mu +{\mu }_{\mathrm{t}})\frac{\partial u}{\partial x})+\frac{\partial }{\partial y}((\mu +{\mu }_{\mathrm{t}})\frac{\partial v}{\partial x})+\frac{\partial }{\partial z}((\mu +{\mu }_{\mathrm{t}})\frac{\partial w}{\partial x})$
v	μ + μt	$S_{\mathrm{v}}=-\frac{\partial p}{\partial y}+\frac{\partial }{\partial x}((\mu +{\mu }_{\mathrm{t}})\frac{\partial u}{\partial y})+\frac{\partial }{\partial y}((\mu +{\mu }_{\mathrm{t}})\frac{\partial v}{\partial y})+\frac{\partial }{\partial z}((\mu +{\mu }_{\mathrm{t}})\frac{\partial w}{\partial y})$
w	μ + μt	$S_{\mathrm{w}}=-\frac{\partial p}{\partial z}+\frac{\partial }{\partial x}((\mu +{\mu }_{\mathrm{t}})\frac{\partial u}{\partial z})+\frac{\partial }{\partial y}((\mu +{\mu }_{\mathrm{t}})\frac{\partial v}{\partial z})+\frac{\partial }{\partial z}((\mu +{\mu }_{\mathrm{t}})\frac{\partial w}{\partial z})$
t	μ/Pr + μt/σt	$S_{\mathrm{t}}$
k	μ + μt/σk	$S_{\mathrm{k}}=\rho G_{\mathrm{k}}-\rho \epsilon$
ε	μ + μt/σε	$S_{\mathrm{\epsilon }}=\frac{\epsilon }{k}(c_1\rho G_{\mathrm{k}}-c_2\rho \epsilon )$

$G_{\mathrm{k}}=\frac{{\mu }_{\mathrm{t}}}{\rho }\left({(\frac{\partial u}{\partial x})}^2+{(\frac{\partial v}{\partial y})}^2+{(\frac{\partial w}{\partial z})}^2+{(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})}^2+{(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})}^2+{(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})}^2\right),$ μ_t=c_μρk²/ε, c_μ=0.09, c₁=1.44, c₂=1.92, σ_t=0.85, σ_ε=1.3, σ_k=1.0.

The supercritical methane-flow in the computational domain is regarded as three-dimensional, incompressible, turbulent and quasi-steady. The boundary conditions for computational model are shown in Figure 4. The flow channels formed by the fins are symmetrical in the longitudinal direction, so the symmetry condition is used for the x-z plane in this study. The flow velocity u_in is assumed to be uniform, the temperature t_in is taken as 195 K, and pressure is specified at the outlet of the channels. On the solid surfaces, no-slip conditions are used, and a constant temperature of 340 K is applied at the top and bottom surfaces of the channels. The temperature distribution in the fins will be determined by solving the conjugated heat transfer problem between supercritical methane and the fins in the computational domain.

Figure 4. Boundary conditions for computational model.

The foregoing governing equations and boundary conditions are solved by Ansys Fluent 19.1. A preprocessor Gambit 2.4 is used to mesh the computational domain for the solver. In order to control the mesh number and improve the grid quality, a hexahedral grid is used for meshing as shown in Figure 5. Due to the stringent requirements for boundary grid density under supercritical flow and heat transfer, five boundary layers are used near the top and bottom walls and fin surfaces, with the thickness of the first layer being 0.01 mm.

Figure 5. Selected grid system.

All variables, including velocity components, pressure and temperature, are averaged to a control volume. The coupling between pressure and velocity is implemented by the semi-implicit method for pressure-linked equations algorithm. The quadratic upstream interpolation for convective kinematics method is used to discretize the convection terms. Palko and Anglart [37] pointed out that the k-ε model failed to calculate the heat transfer to supercritical fluid in deteriorated region, but the shear stress transport (SST) k-ω turbulence model predicted the heat transfer coefficient in a very good agreement with the experimental results. Therefore, the SST k-ω turbulence model is adopted to simulate the fluid flow. The SST k-ω turbulence model [38] is a hybrid model combining the Wilcox k-ω and the k-ε models. A blending function, F1, activates the Wilcox model near the wall and the k-ε model in the free stream. This ensures that the appropriate model is utilized throughout the flow field: The k-ω model is well-suited for simulating flow in the viscous sublayer, and the k-ε model is ideal for predicting flow behavior in regions away from the wall. The residuals are set to be less than 10⁻⁵ and 10⁻⁸ for continuity and energy equations, respectively, to ensure convergence of the computations.

In order to validate the solution independency of the grid number, four different grid systems for different cases are investigated, and the predicted average Nusselt numbers and f-factors are shown in Table 6. Compared with the result of the finest grid systems for different cases, the deviations of the average Nusselt number and f-factor of the second finest grid systems are less than 2.0%, therefore, the reference grid systems for Case A, Case B and Case C–F are about 554,600 cells, 833,600 cells and 1,132,000 cells, respectively.

Table 6. Variation of grid independence (u_in=8.0 m/s).

Configuration	Grid	$\overline{Nu}$	Relative error	f	Relative error
Case A	277300	184.3	5.40%	0.019156	−6.80%
	415950	180.8	3.40%	0.019773	−3.80%
	554610	177.16	1.30%	0.02025	−1.50%
	693250	174.86	–	0.020554	–
Case B	415930	259.19	6.70%	0.018982	−7.60%
	623900	252.63	4.00%	0.01966	−4.30%
	833620	246.61	1.50%	0.02024	−1.50%
	1031870	242.91	–	0.020544	–
Case C	452830	337.98	7.70%	0.037283	−8.60%
	735860	327.62	4.40%	0.038548	−5.50%
	1132170	319.24	1.70%	0.04007	−1.80%
	1471530	313.81	–	0.040791	–

Data reduction and code validation

The Reynolds number (Re) is defined as (2) $$Re=\frac{\rho u_{\text{max}}{D}_{\mathrm{h}}}{\mu }=\frac{m{D}_{\mathrm{h}}}{\mu A_{\min }}$$(2) where u_max is the maximum velocity of methane evaluated at the minimum cross section, m is the mass flow rate, and A_min is the minimum cross section in the channels. ρ is the density of methane, and μ is the dynamic viscosity of methane, which are evaluated at the mean temperature of methane. D_h is the hydraulic diameter and is defined by

(3) $$D_{\mathrm{h}}=\frac{4A_{\min }}{P}$$(3) where P is the perimeter.

The local Nusselt number (Nu_L) is defined as (4) $$N{u}_{\mathrm{L}}=\frac{h_{\mathrm{L}}{D}_{\mathrm{h}}}{{\lambda }_{\mathrm{L}}}$$(4)

The mean Nusselt number ($\overline{Nu}$) is defined as [22] (5) $$\overline{Nu}=\frac{1}{n}\sum _{i=1}^{n}N{u}_{\mathrm{L}}=\frac{1}{n}\sum _{i=1}^n\frac{m{c}_{\mathrm{p},\mathrm{i}}\left(t_{\mathrm{i}+1}-t_{\mathrm{i}}\right)D_{\mathrm{h}}}{\left(A_{\mathrm{f},\mathrm{i}}{\eta }_{\mathrm{f},\mathrm{i}}+A_{\mathrm{w},\mathrm{i}}\right)\Delta T_{\mathrm{i}}{\lambda }_{\mathrm{i}}}=\frac{1}{n}\sum _{i=1}^n\frac{Q_{\mathrm{i}}{D}_{\mathrm{h}}}{{\eta }_{\mathrm{o},\mathrm{i}}{A}_{\mathrm{o},\mathrm{i}}\Delta T_{\mathrm{i}}{\lambda }_{\mathrm{i}}}$$(5) in which, the local logarithmic mean temperature difference ΔT_i is defined as [39] (6) $$\Delta T_{\mathrm{i}}=\frac{t_{\mathrm{i}+1}-t_{\mathrm{i}}}{\hspace{0.17em}\text{ln}\hspace{0.17em}\left(\raisebox{1ex}{$\left(t_{\mathrm{w}}-t_{\mathrm{i}}\right)$}\!\left/ \!\raisebox{-1ex}{$\left(t_{\mathrm{w}}-t_{\mathrm{i}+1}\right)$}\right.\right)}$$(6) where h_L is the local heat transfer coefficient, n is the total cross section number. Q_i is the local heat transfer rate, t_i is the temperature of methane in the i cross section, t_w is the temperature of the top and bottom surfaces of the local passage, while A_w,i is the corresponding area. A_f,i is the fin area, and A_o,i is the total heat transfer area of the local passage, A_o,i= A_f,i+ A_w,i. η_f,i is the fin efficiency, and η_o,i is the surface efficiency. λ_i is the thermal conductivity of methane and c_p,i is the specific heat of methane, which are evaluated at the mean temperature of methane in the local passage. The fin efficiency η_f,i and the surface efficiency η_o,i are described as (7) $${\eta }_{\mathrm{f},\mathrm{i}}=\frac{Q_{\mathrm{i}}}{Q_{\text{ideal},\mathrm{i}}}$$(7) (8) $${\eta }_{\mathrm{o},\mathrm{i}}=1-\frac{A_{\mathrm{f},\mathrm{i}}}{A_{\mathrm{o},\mathrm{i}}}\left(1-{\eta }_{\mathrm{f},\mathrm{i}}\right)$$(8) where Q_ideal is the ideal heat transfer rate at the condition of the whole fin temperature equal to t_w.

The Colburn factor (j-factor) is defined as (9) $$j=\frac{\overline{Nu}}{\mathit{ReP}{r}^{1/3}}$$(9)

The Darcy friction factor (f-factor) can reflect the pressure drop of the channel and is given as (10) $$f=\frac{2\Delta p_{\mathrm{f}}}{\rho {\left(u_{\text{max}}\right)}^2}\frac{D_{\mathrm{h}}}{L_2}$$(10) where Δp_f is the frictional pressure drop and L₂ is the length of the channels. The frictional pressure drop Δp_f can be obtained by (11) $$\Delta p_{\mathrm{f}}=\Delta p-\Delta p_{\mathrm{a}}$$(11) (12) $$\Delta p_{\mathrm{a}}={\rho }_{\text{out}}{\left(u_{\text{out}}\right)}^2-{\rho }_{\text{in}}{\left(u_{\text{in}}\right)}^2$$(12) where Δp is the total pressure drop, and Δp_a is the acceleration pressure drop. ρ_in, u_in and ρ_out, u_out are the density and velocity evaluated at the inlet and outlet of the channels, respectively.

In order to validate the reliability of the present numerical method, a model set up from the experiments of Dang and Hihara [40] is simulated and compared with the experimental correlations near the pseudo-critical temperature. The thermal and hydraulic correlations used for comparison are described as follows. (13) $$\overline{Nu}=\frac{\left(f_{\mathrm{f}}/8\right)\left(R{e}_{\mathrm{b}}-1000\right)Pr}{1.07+12.7\sqrt{f_{\mathrm{f}}/8}\left(P{r}^{2/3}-1\right)}$$(13) (14) $$Pr=\{\begin{array}{c}{c}_{\mathrm{p},\mathrm{b}}{\mu }_{\mathrm{b}}/{\lambda }_{\mathrm{b}},\text{for}{c}_{\mathrm{p},\mathrm{b}}\ge {\overline{c}}_{\mathrm{p}}\\ {\overline{c}}_{\mathrm{p}}{\mu }_{\mathrm{b}}/{\lambda }_{\mathrm{b}},\text{for}{c}_{\mathrm{p},\mathrm{b}}<{\overline{c}}_{\mathrm{p}}\text{and}{\mu }_{\mathrm{b}}/{\lambda }_{\mathrm{b}}\ge {\mu }_{\mathrm{f}}/{\lambda }_{\mathrm{f}}\\ {\overline{c}}_{\mathrm{p}}{\mu }_{\mathrm{f}}/{\lambda }_{\mathrm{f}},\text{for}{c}_{\mathrm{p},\mathrm{b}}<{\overline{c}}_{\mathrm{p}}\text{and}{\mu }_{\mathrm{b}}/{\lambda }_{\mathrm{b}}<{\mu }_{\mathrm{f}}/{\lambda }_{\mathrm{f}}\end{array}$$(14) (15) $${\overline{c}}_{\mathrm{p}}=\frac{h_{\mathrm{w}}^{\ast }-h_{\mathrm{b}}^{\ast }}{t_{\mathrm{w}}-t_{\mathrm{b}}}$$(15) (16) $$R{e}_{\mathrm{b}}=\frac{Gd}{{\mu }_{\mathrm{b}}}$$(16) (17) $$f_{\mathrm{f}}={\left(1.82{\hspace{0.17em}\text{log}\hspace{0.17em}}_{10}\left(R{e}_{\mathrm{f}}\right)-1.64\right)}^{-2}$$(17) (18) $$R{e}_{\mathrm{f}}=\frac{Gd}{{\mu }_{\mathrm{f}}}$$(18) where subscript b represents the mean bulk temperature, and subscript f represents the film temperature. t_f can be obtained by (19) $$t_{\mathrm{f}}=\frac{t_{\mathrm{b}}+t_{\mathrm{w}}}{2}$$(19)

In above correlations, the range of mass flux G is from 200 to 1200 kg/(m²·s), the range of inner diameter d of tube is from 1 to 6 mm, and Prandtl number Pr ranges from 1.3 to 12.0.

Figure 6 shows the comparisons between the simulated and experimental results. The simulated results of the SST k-ω turbulence model agree well with the experimental results. Compared with the experimental results, the maximum and mean relative deviations of the average Nusselt number are 13.5 and 6.7%, respectively, and the maximum and mean relative deviations of the f-factor are 7.8% and 6.2%, respectively. Such good agreements between the simulated and experimental results secure the reliability of the present numerical method.

Figure 6. $\overline{Nu}$ and f-factor comparisons between simulated and experimental results.

Numerical results and discussions

Effect of pressure on the thermal-hydraulic performance

The working conditions are summarized in Table 7 in this study. The pseudo-critical temperature of methane at 7, 8, 9, and 10 MPa is about 205.2, 210.1, 214.3, and 218.1 K, respectively. From Figure 1, it can be seen that the temperature of methane along the flow direction in PCHE channels for different cases goes through the pseudo-critical temperature at different operating pressures.

Table 7. Working conditions for different cases.

Parameter	Symbol/unit	Value
inlet temperature	tin / K	195
temperature difference	Δt / K	>30
inlet pressure	pin / MPa	7–10
pressure drop	Δpf / MPa	<0.1
mass flux	G / kg/m^2⋅s	132.32–2117.12
wall temperature	tw / K	340

In a PCHE channel, the temperature of supercritical methane will gradually increase in the heating process and the physical properties will undergo dramatic changes. The comparisons of specific heat and density at the central cross section in the z-direction for Case B at different pressures are shown in Figures 7 and 8. From Figure 7, it can be seen that the specific heat increases first and then decreases with the increase of temperature, and the specific heat decreases with the increase of pressure at the same location in the channel. The specific heat of supercritical methane is like a pulse function near the pseudo-critical point in Figure 1, and large specific heat can significantly improve the heat transfer performance. As can be seen from Figure 8, the density of supercritical methane increases with the increase of pressure, so the velocity at the same mass flux is higher at lower operating pressure, then the turbulence intensity of the main flow is larger inducing the heat transfer enhancement.

Figure 7. Comparisons of c_p for Case B at different pressures.

Figure 8. Comparisons of ρ for Case B at different pressures.

The local thermal and hydraulic performance of Case B at different operating pressures along the main flow direction is described in Figure 9. As shown in Figure 9, the local Nusselt number Nu_L in the first row is the highest due to the inlet effect, and then the local Nusselt number Nu_L decreases along the main flow direction and reduces to the minimum in the second row. After that, it increases and achieves approximately fully developed condition since the 6th row. The curves of specific heat at different operating pressures show a parabolic distribution along the main flow direction, and specific heat reaches its peak value at the pseudo-critical point. From Figure 9, it can be seen that the local heat transfer performance is reduced as pressure increases, but it is worth to point out that the local heat transfer performance is enhanced near the pseudo-critical point at different operating pressures.

Figure 9. Local Nusselt number Nu_L and c_p along the main flow direction at u_in=2.0 m/s.

The average thermal and hydraulic performance of Case B is compared at different pressures as shown in Figure 10, where the mass flux and Reynolds number in the channel is taken as the abscissa. In Figure 10(a), it is found that the average Nusselt number $\overline{Nu}$ increases with the increase of mass flux, and the heat transfer performance is reduced as pressure increases. This can be explained as follows. At the same mass flux, the methane density in the channel is smaller at lower pressure as can be seen in Figure 8, so the methane velocity is higher at lower pressure, and the higher velocity can improve the heat transfer performance. On the other hand, the specific heat is larger at lower pressure in Figure 7, and the larger specific heat can greatly improve the capability of supercritical methane for heat transfer enhancement. As a consequence of the above two main reasons, the heat transfer performance is increased as the pressure decreases. From Figure 10(b), it can be seen that the Reynolds number is different at the same mass flow rate, and under the same mass flow rate, the Reynolds number at 7 MPa is the largest among these four different operating pressures. This is because that the terms m, D_h, and A_min are the same in Eq. (2)(2) $$Re=\frac{\rho u_{\text{max}}{D}_{\mathrm{h}}}{\mu }=\frac{m{D}_{\mathrm{h}}}{\mu A_{\min }}$$(2) , but the value of μ at 7 MPa is the largest as can be seen from Figure 1(d). Figure 10 also shows the variation of the f-factor with the change of mass flux at different operating pressures. With the increase of mass flux, the f-factor decreases. The reason is that the methane velocity is decreased with the increase of pressure, the friction performance is reduced as the operating pressure increases from 7.0 to 10.0 MPa.

Figure 10. $\overline{Nu}$ and f-factor comparisons for Case B at different pressures.

Thermal-hydraulic performance comparison for different PCHE configurations

The inlet pressure for Cases A–F is set at 8 MPa, and the inlet velocities are set at 0.5, 0.7, 1.0, 2.0, 4.0, 6.0, and 8.0 m/s, respectively. The heat transfer characteristics of different PCHE configurations are compared in Figure 11, where the inlet velocity and Reynolds number in channels is taken as the abscissa. As can be seen from Figure 11, the heat transfer performance of different PCHEs increases with the increase of inlet velocity and Reynolds number. At the same inlet velocity, Case C provides the best effectiveness of the heat transfer enhancement, the next is Case D, then Case F, Case E, Case B in descending order, and the last one is Case A. The comparative results of all delta VG and airfoil configurations can indicate that delta VGs arranged upstream the airfoils can provide better heat transfer performance than delta VGs arranged downstream the airfoils, and the heat transfer performance of the common-flow-up configuration is not always better than that of the common-flow-down configuration. The Reynolds number is different at the same inlet velocity, and the value of Case C is the largest. The reason can be explained as follows. Case C provides the best heat transfer performance, and the mean temperature in the PCHE channel of Case C is the largest. From Figure 1(d), it can be seen that the value of μ at larger mean temperature is smaller in the studied temperature range, so the Reynolds number of Case C is the largest according to Eq. (2)(2) $$Re=\frac{\rho u_{\text{max}}{D}_{\mathrm{h}}}{\mu }=\frac{m{D}_{\mathrm{h}}}{\mu A_{\min }}$$(2) .

Figure 11. Heat transfer performance comparisons for different PCHEs.

Figure 12 shows the streamlines at the central cross section in the z-direction of Case B at u_in=2.0 m/s. From Figure 12, it can be observed that the flow direction of supercritical methane upstream the airfoil is almost parallel to the x-direction, then the flow direction is changed and consistent with the direction of the airfoil surface. Figure 13 shows the streamlines in the 7th row of Cases C–F along the main flow direction at u_in=2.0 m/s. As can be seen from Figure 13, the flow direction of supercritical methane upstream the VG is changed, and LVs are generated downstream the VG when supercritical methane flows over the VG, and it can destroy the boundary layer effectively. It can be also seen that Cases C–F can provide higher velocity along the airfoil fin than Case B, especially Cases C and D. As is well known, the higher velocity can greatly improve the capability of supercritical methane for heat transfer enhancement. Therefore, VGs arranged upstream the airfoils, that is, Case C and Case D, can effectively enhance heat transfer. For Case E, VGs arrangement is almost consistent with the flow direction in Figure 13(c), so the methane flow behind the VG cannot be effectively disturbed, and the heat transfer enhancement is worse than Case F. That is to say, the arrangement of VGs should be properly designed in order to further enhance heat transfer.

Figure 12. Streamlines at the central cross section in the z-direction of Case B at u_in = 2.0 m/s.

Figure 13. Streamlines in the 7th row of Cases C–F at u_in = 2.0 m/s.

To avoid the effect of the entrance or the exit, the 7th airfoil fin is picked out and its local performance is analyzed. There are 16 equidistant cross-sectional planes from x = 48 mm plane to x = 56 mm plane. Figure 14 presents the local Nusselt number Nu_L at the cross-sectional planes along the flow direction for Cases B–F at u_in=2.0 m/s. It can be seen that the local heat transfer performance of channels with airfoil and VG fins is always higher than that of the airfoil channel behind the trailing edge of VGs, and the main difference of the local Nusselt number Nu_L among Cases C–F can be attributed to the VG layout. As mentioned above, the intensity of LVs generated by delta VGs arranged upstream the airfoils is stronger than that of LVs generated by delta VGs arranged downstream the airfoils, therefore, the local Nusselt number Nu_L of Cases C and D is higher than that of Cases E and F.

Figure 14. Local Nusselt number Nu_L in the 7th row at u_in=2.0 m/s.

Figure 15 presents the velocity distributions at the different cross sections in the y-z plane of Cases C–F at u_in=2.0 m/s. From Figure 15, it can be observed that LVs are generated at the leading edges of Cases C–F when supercritical methane flows over the VGs. The structure of NACA 23024 airfoil is not symmetrical, so the left LVs and right LVs are also asymmetrical. LVs occur not only in the region where the VGs are mounted, but also survive in a long distance downstream. The interaction between adjacent vortices can be divided into two types, that is, counter-rotating and co-rotating. Counter-rotating vortices are identified by the direction of the secondary flow between the pairs which can be directed either toward or away from the wall. For common-flow-down configuration (Case D and Case F), the transverse distance between the leading edges of the VG pairs is less than the transverse distance between the trailing edges, while, for the common-flow-up configuration (Case C and Case E), the transverse distance between the leading edges is larger than that of the trailing edges. The counter-rotating pairs of vortices are generated by Case D and Case F. The left vortices rotate in a clockwise direction, while the right vortices rotate in a counterclockwise direction. These counter-rotating vortices create a down-wash flow toward the lower channel wall between the vortices and the up-wash flow away from the lower wall in the region outside of the vortices. The counter-rotating pairs of vortices are also generated by Case C and Case E but the left vortices rotate in a counterclockwise direction, while the right vortices rotate in a clockwise direction. These vortices are responsible for the up-wash flow between the vortices and the down-wash flow in the outer region of the vortices. These two vortices both induce strong swirling motion and destroy the boundary layer, therefore, they may greatly enhance heat transfer. When LVs, generated by delta VGs arranged upstream the airfoils, flow through the airfoils, the intensity of LVs is stronger than that of LVs generated by delta VGs arranged downstream the airfoils. The LVs generated by Case C and Case D have greater influence downstream, therefore, delta VGs arrangement upstream the airfoils is a good choice for heat transfer enhancement.

Figure 15. Velocity distributions at the different cross-sections at u_in=2.0 m/s.

Figure 16 illustrates the variations of friction factor for all cases above. The f-factor decreases with the increase of inlet velocity. The comparative results indicate that the f-factors of Case C and Case E are smaller than those of Case D and Case F at the same inlet velocity, respectively. The results show that the pressure drop performance of the common-flow-up configuration is better than that of the common-flow-down configuration. The reason can be explained as follows. The flow resistance of supercritical methane from the inlet to outlet of the channel with VGs and airfoils is composed by the acceleration resistance of supercritical methane, the frictional resistance of the channel surfaces and the local resistance of the VGs and airfoils, and the local resistance of the VGs and airfoils plays a dominant part of the total pressure drop. The skin friction perturbation of the common-flow-down cases is much larger than that of the common-flow-up cases when the fluid flows through the VG pairs. This is so because that there is a stronger interaction effect on the wall. The spanwise component of the skin friction of the common-flow-down cases is much larger, and the spanwise skin friction of the common-flow-up cases will decrease even more as the vortices continue to lift each other away from the wall. As a consequence, the local pressure drop along common-flow-down VG pairs is larger, that is, the hydraulic performance of Case D and Case F is worse than that of Case C and Case E, respectively.

Figure 16. Friction factor performance comparisons for different PCHE configurations.

Overall heat transfer performance comparison

As mentioned above, Case C provides the best heat transfer performance, and Case A provides the smallest pressure drop. Attention is now turned to the overall heat transfer performance of different PCHE configurations at the same mass flow rate. In the present study, the volume goodness criterion is adopted.

The volume goodness criterion was suggested by Kays and London [41]. The volume goodness factors are given as follows (20) $${\eta }_{\mathrm{o}}{h}_{\text{std}}\beta =\frac{c_{\mathrm{p}}\mu }{P{r}^{2/3}}{\eta }_{\mathrm{o}}\frac{4\sigma }{D_{\mathrm{h}}{}^2}jRe$$(20) (21) $$E_{\text{std}}\beta =\frac{{\mu }^3}{2{\rho }^2}\frac{4\sigma }{D_{\mathrm{h}}{}^4}fR{e}^3$$(21) where σ is the ratio of the minimum flow area A_min to the frontal area A_fr. η_oh_stdβ represents the heat transfer power per unit temperature difference and per unit core volume, and E_stdβ represents the friction power expenditure per unit core volume. The subscripts std is any preferred standard condition for comparison. A higher plot on η_oh_stdβ vs. E_stdβ characterizes that a surface is better from the viewpoint of heat exchanger volume required. Geometries having large values of η_oh_stdβ will require less core volume for a given fluid-side performance. As shown in Figure 17, Case C has the best volume goodness factor, the next is Case E, then Case F, Case D, Case B in descending order, and the last one is Case A.

Figure 17. Volume goodness comparisons for different PCHE configurations.

Conclusions

This article compares the effects of six different PCHE configurations on the heat transfer and flow friction characteristics for supercritical LNG by three dimensional numerical simulations. Main conclusions can be summarized as follows:

1. At the same inlet velocity, the effects of different heat transfer enhancement models are concluded such that Case C > Case D > Case F > Case E > Case B > Case A. The comparative results of all delta VG and airfoil configurations can indicate that delta VGs arranged upstream the airfoils can provide better heat transfer performance than delta VGs arranged downstream the airfoils. Heat transfer performance of common-flow-up configuration is not always better than that of common-flow-down configuration. This means that the arrangement of VGs should be properly designed in order to further enhance heat transfer.

2. Moreover, the comparative results indicate that at the same inlet velocity, the f-factors of Case C and Case E are smaller than those of Case D and Case F, respectively. The results of the numerical simulations depict that the hydraulic performance of the common-flow-up configuration is better than that of the common-flow-down configuration.

3. The volume goodness criterion is adopted to evaluate the overall heat transfer performance of the different PCHE configurations at the same mass flow rate. The results show that Case C has the best volume goodness factor, the next is Case E, then Case F, Case D, Case B in descending order, and the last one is Case A.

Nomenclature
A	=	area, m2
Af	=	fin area, m2
Afr	=	frontal area, m2
Amin	=	minimum flow area, m2
Ao	=	total surface area, m2
Aw	=	wall surface area, m2
c1, c2, cμ	=	turbulence model constant
cp	=	specific heat at constant pressure, J/(kg⋅K)
${\overline{c}}_{\mathrm{p}}$	=	${\overline{c}}_{\mathrm{p}}=(h_{\mathrm{w}}^{*}-h_{\mathrm{b}}^{*})/(t_{\mathrm{w}}-t_{\mathrm{b}})$
d	=	inner diameter, mm
Dh	=	hydraulic diameter, mm
Estdβ	=	friction power expenditure per unit core volume, see Eq. (21)(21) $$E_{\text{std}}\beta =\frac{{\mu }^3}{2{\rho }^2}\frac{4\sigma }{D_{\mathrm{h}}{}^4}fR{e}^3$$(21)
f	=	Darcy friction factor
G	=	mass flux, kg/(m2⋅s)
Gk	=	generation of turbulence kinetic energy, kg/(m⋅s3)
h	=	heat transfer coefficient, W/(m2⋅K)
h*	=	enthalpy, J/kg
H	=	channel height, mm
Hf	=	airfoil height, mm
HVG	=	VG height, mm
j	=	Colburn factor
k	=	kinetic energy, m2/s2
L1, L2, L3	=	channel length, mm
Lf	=	airfoil chord length, mm
LNG	=	liquefied natural gas
LNG- FPSO	=	liquefied natural gas - floating production storage and offloading
LNG- FSRU	=	liquefied natural gas - floating storage and regasification units
LV	=	longitudinal vortex
LVG	=	VG length, mm
m	=	mass flow rate, kg/s
n	=	total cross section number
N	=	number of fin rows
NACA	=	National Advisory Committee for Aeronautics
NG	=	natural gas
Nu	=	Nusselt number
$\overline{Nu}$	=	mean Nusselt number
p	=	pressure, Pa
P	=	perimeter, mm
PCHE	=	printed circuit heat exchanger
P1	=	longitudinal pitch, mm
Pt	=	transverse pitch, mm
Pr	=	Prandtl number
Q	=	heat transfer rate, W
Re	=	Reynolds number
SST	=	shear stress transport
Sϕ	=	source term, W/m3
t	=	temperature, K
u, v, w	=	x, y, z velocity components, m/s
UDF	=	user defined function
VG	=	vortex generator
W	=	channel width, mm
Wf	=	maximum airfoil thickness, mm
WVG	=	VG thickness, mm
x, y, z	=	Cartesian co-ordinates

Greek symbols
α	=	VG attack angle, °
Γϕ	=	diffusion coefficient, m2/s
Δp	=	total pressure drop, Pa
Δpa	=	acceleration pressure drop, Pa
Δpf	=	frictional pressure drop, Pa
Δt	=	temperature difference, K
ΔT	=	logarithmic mean temperature difference, K
ε	=	dissipation rate, 1/s
ηf	=	fin efficiency
ηo	=	surface efficiency
ηohstdβ	=	heat transfer power per unit temperature difference and per unit core volume, see Eq. (20)(20) $${\eta }_{\mathrm{o}}{h}_{\text{std}}\beta =\frac{c_{\mathrm{p}}\mu }{P{r}^{2/3}}{\eta }_{\mathrm{o}}\frac{4\sigma }{D_{\mathrm{h}}{}^2}jRe$$(20)
λ	=	thermal conductivity, W/(m⋅K)
μ	=	dynamic viscosity, kg/(m⋅s)
μt	=	turbulent eddy-viscosity, kg/(m⋅s)
ρ	=	density, kg/m3
σ	=	ratio of the minimum flow area Amin to the frontal area Afr
σk, σt, σε	=	turbulence model constant
ϕ	=	dependent variable
ω	=	specific dissipation rate, 1/s

Subscripts
b	=	mean value
f	=	fin, frictional, film temperature
fr	=	frontal area
i	=	i cross section
ideal	=	ideal condition
in	=	inlet
L	=	local
max	=	maximum value
min	=	minimum value
out	=	outlet
std	=	standard condition
w	=	wall

Additional information

Funding

This work is supported by National Nature Science Foundation of China (Grant No. 51774237) and The Youth Scientific Research and Innovation Team of Xi’an Shiyou University (No. 2019QNKYCXTD10).

Notes on contributors

Ling-Hong Tang

Linghong Tang is a visiting scholar at the Department of Energy Sciences of Lund University. He is an associate professor in the School of Mechanical Engineering, Xi’an Shiyou University. He received his Ph.D. degree in engineering thermophysics from Xi’an Jiaotong University in 2009, and he finished his postdoctoral research from Harbin Engineering University in 2014. His main research interests are numerical heat transfer and compact heat exchangers.

Jie Pan

Jie Pan is an associate professor at the Department of Oil and Gas Storage and Transportation Engineering of Xi’an Shiyou University. His major research fields are multiphase flow in petroleum engineering, LNG heat exchanger and LNG cold energy utilization.

Bengt Sundén

Bengt Sundén received his M.Sc. in Mechanical Engineering 1973, Ph.D. in Applied Thermodynamics and Fluid Mechanics 1979, became Docent in Applied Thermodynamics and Fluid Mechanics 1980, all at Chalmers University of Technology, Göteborg, Sweden. He is a professor of Heat Transfer, Lund University, and he serves as Professor Emeritus and Senior Professor in Heat Transfer since 2016-10-01. He is a fellow of the ASME. The research activities cover compact heat exchangers, heat transfer enhancement, gas turbine heat transfer, combustion-related heat transfer including thermal radiation, CFD-methods, liquid crystal thermography, condensation and evaporation on micro-nanostructured surfaces, nanofluids, impinging jets, aerospace heat transfer, computational modeling and analysis of multiphysics and multiscale phenomena for fuel cells (SOFC, SOEC, PEMFC), thermal management of batteries.