*2.1. Geometry and Boundary Conditions*

The flow and heat transfer performance of CO2 in horizontal zigzag channels is investigated in this paper. Figure 2 shows a schematic diagram of the computational domain with the boundary condition adopted in this research. The study will consider the effects of different geometric factors, including hydraulic diameter, bending angle, and pitch length on the internal flow and heat transfer characteristics of the channel, as shown in Table 1. To check the independent effect of the geometric parameter, the comparative study only changes one geometric parameter at a time, and the rest remains unchanged. For example, for investigating the effect of bend angle *θ* on the channel flow and heat transfer performance, four types of bend angle, 100◦, 115◦, 140◦, and 180◦ (straight), are considered for the geometric model, while *Lp* is fixed to 7.75 mm and *dh* is fixed to 2 mm.

As shown in Table 2, the inlet temperature was changed from 280 K to 350 K to ensure the bulk temperature *Tb* covers the pseudocritical point of CO2 for this analysis. The outlet pressure of CO2 varies between 7.5 MPa and 9 MPa to keep its condition near the pseudocritical point. Constant wall heat flux (12 kW/m2 for heating case and −12 kW/m2 for cooling case) and mass flux 200 kg/m2-s were adopted in this numerical simulation. The thermal properties of CO2 were obtained from NIST Database (REFPROP V9.1).

**Figure 2.** Physical model of zigzag channel.

**Table 1.** Geometric parameters.


**Table 2.** Boundary conditions in detail.


#### *2.2. Governing Equations*

The governing equations regarding the continuity, momentum, and energy are expressed as Equations (1)–(3).

(1) Continuity equation:

$$\frac{\partial}{\partial \mathbf{x}\_i}(\rho u\_i) = 0 \tag{1}$$

(2) Momentum equation [25]:

$$\frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} (\rho u\_i u\_{\dot{j}}) = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial \tau\_{\dot{i}j}}{\partial \mathbf{x}\_{\dot{j}}} + \rho g\_i \tag{2}$$

(3) Energy equation:

$$\frac{\partial}{\partial \mathbf{x}\_i} (\mu\_i (\rho E + p)) = \frac{\partial}{\partial \mathbf{x}\_i} \left( k\_{eff} \frac{\partial T}{\partial \mathbf{x}\_i} \right) \tag{3}$$

where *ρgi* is the gravitational body force, *ui* is overall velocity vector, *T* is the temperature, *ρ* is the density, *ui* is the velocity vector, *E* is the total energy, *p* is the static pressure, *τij* is the stress tensor, and *keff* is the effective conductivity (*keff* = *k* + *kt*, *kt* is the turbulent thermal conductivity).

(4) Direct numerical simulation (DNS) of Navier–Stokes equations is the most accurate method for turbulence; however, it is not feasible in most situations to resolve the wide range of scales in time and space as the CPU requirements by far exceed the existing capacity. For this reason, averaging procedures such as the Reynolds method have to be applied to the Navier–Stokes equations to filter out the turbulent spectrum [26]. However, the averaging process introduces additional unknown terms into the transport equations that need to be provided by suitable turbulence models. The SST *k*–*ω* turbulence model is used in this paper and is described as follows.

The transport equations of the *k*–*ω* model are expressed as Equations (4) and (5):

$$\frac{\partial}{\partial \mathbf{x}\_i} (\rho k u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left( \Gamma\_k \frac{\partial k}{\partial \mathbf{x}\_j} \right) + G\_k - \rho \beta^\* k \omega \tag{4}$$

$$\frac{\partial}{\partial \mathbf{x}\_i} (\rho \omega u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left( \Gamma\_w \frac{\partial \omega}{\partial \mathbf{x}\_j} \right) + \mathbf{G}\_{\omega} - \rho \beta \omega^2 \tag{5}$$

where *Gk* is the generation of turbulence kinetic energy *k* due to mean velocity gradients, *G<sup>ω</sup>* is the generation of *ω*, and Γ*<sup>k</sup>* and Γ*<sup>w</sup>* represent the effective diffusivity of *k* and *ω* calculated by Equations (6) and (7).

$$
\Gamma\_k = \mu + \mu\_t \left(\frac{F\_1}{\sigma\_{k1}} + \frac{1 - F\_1}{\sigma\_{k2}}\right) \tag{6}
$$

$$
\Gamma\_{\omega} = \mu + \mu\_t \left( \frac{F\_1}{\sigma\_{\omega 1}} + \frac{1 - F\_1}{\sigma\_{\omega 2}} \right) \tag{7}
$$

where *σ<sup>k</sup>* and *σω* are turbulent Prandtl numbers for *k* and *ω*, respectively, *μ* is the dynamic viscosity of the fluid, and *F*<sup>1</sup> is calculated as Equations (8)–(10).

$$F\_1 = \tanh\left(\phi\_1^4\right) \tag{8}$$

$$\phi\_1 = \min\left[ \max\left( \frac{\sqrt{k}}{0.09\omega y'}, \frac{500\mu}{\rho y^2 \omega} \right), \frac{4\rho k}{\sigma\_{\omega 2} D\_\omega y^2} \right] \tag{9}$$

$$D\_{\omega} = 2(1 - F\_1)\rho \sigma\_{\omega 2} \frac{1}{\omega} \frac{\partial k}{\partial x\_j} \frac{\partial \omega}{\partial x\_j} \tag{10}$$

The turbulent viscosity *μ<sup>t</sup>* of the SST *kω* model is calculated using Equations (11)–(15).

$$\mu\_t = \frac{\rho k}{\omega} \frac{1}{\max\left(\frac{1}{a^\pi}, \frac{\sqrt{2S\_{ij}}F\_2}{a\_1\omega}\right)}\tag{11}$$

$$S\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{12}$$

$$F\_2 = \tanh\left(\phi\_2^2\right) \tag{13}$$

$$\phi\_2 = \max\left(2\frac{\sqrt{k}}{0.09\omega y'}\frac{500\mu}{\rho y^2 \omega}\right) \tag{14}$$

$$a^\* = a^\*\_{\infty} \left( \frac{a^\*\_0 + \frac{\rho k}{6\mu \omega}}{1 + \frac{\rho k}{6\mu \omega}} \right) \tag{15}$$

where *μ<sup>t</sup>* is the turbulent viscosity, *y* is the wall distance, and the constants used in the SST model are as follows: *a*∗ <sup>∞</sup> = 1, *a*1= 0.31, *σk*<sup>1</sup> = 1.176, *σk*<sup>2</sup> = 1, *σω*<sup>1</sup> = 2, *σω*<sup>2</sup> = 1.168, *a*<sup>∗</sup> <sup>0</sup> = 0.024.

The second-order upwind scheme of Equation (16) is used to discretize the convective term of the above governing Equation [27].

$$
\varphi\_{f,su} = \varphi + \nabla \varphi \cdot \vec{r} \tag{16}
$$

where *<sup>ϕ</sup>* and <sup>∇</sup>*<sup>ϕ</sup>* are the cell-centered value and its gradient in the upstream cell and <sup>→</sup> *r* is the displacement vector from upstream cell centroid to the face centroid.

*2.3. Data Reduction*

The hydraulic diameter, *dh*, was defined as Equation (17):

$$d\_h = \frac{4A}{\mathbb{C}\_{\text{net}}}\tag{17}$$

where *A* is the semicircular cross section area of the channel, and *Cwet* is the wet circumference of the cross section.

The average heat convection coefficient, *h*, was defined as Equation (18):

$$h = \frac{Q\_w}{T\_w - T\_b} \tag{18}$$

where *Qw* is the heat flux of the wall, *Tw* is the average wall temperature, and *Tb* is the fluid bulk temperature.

The channel total pressure drop, Δ*P*, was defined as Equation (19):

$$
\Delta P = P\_{\rm in} - P\_{\rm out} \tag{19}
$$

where *Pin* and *Pout* are area-weighted average pressure at the inlet and outlet of the channel, respectively.

#### *2.4. Grid Independence and Model Validation*

Meshes in this study are generated using ICEM CFD 2019, and the size of the first layer adjacent to the wall is less than 2 × <sup>10</sup>−<sup>6</sup> <sup>m</sup> to keep the wall y<sup>+</sup> < 1. Structured hexahedral cells are used in the whole computational domain, and the mesh is locally densified at the bend of the flow direction, as illustrated in Figure 3.

**Figure 3.** Mesh structure of the computational domain.

Mesh independence analysis was performed with five different mesh sizes of 91,656, 152,988, 233,508, 314,028, and 394,548. Percentage changes in physical variables *h* and Δ*P* were chosen as the basis for independence judgment. As the comparison result listed in Table 3, the relative error percentage of *h* changes from 11.01% to 1.17%, and the relative error percentage of Δ*P* changes from 19.78% to 0.59% with the mesh refinement. Consequently, the mesh size of 314,028 elements is considered sufficient, and 300,000 is used as the baseline for all the rest of the simulations.


**Table 3.** Mesh independence analysis.
