CFD Modelling of Pressure Drop in Double-Pipe Heat Exchanger with Turbulators Using a Porous Media Model

The main governing equations for CFD simulations are conservation equations for mass (Eq. (1)) and momentum (Eq. (2)):

$-\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \vec{v})=0$        (1)

$\frac{\partial}{\partial t}(\rho \vec{v})+\nabla \cdot(\rho \vec{v} \vec{v})=-\nabla p+\nabla \cdot(\bar{\bar{\tau}})+\rho \vec{g}+\vec{F}$        (2)  

$\bar{\bar{\tau}}=\mu\left[\left(\nabla \vec{v}+\nabla \vec{v}^T\right)-\frac{2}{3} \nabla \cdot \nabla \vec{v} I\right]$        (3)

2.1 Turbulence models

The exact solution of the governing Navier-Stokes equations for turbulent flows is possible only for simple single-phase cases; it requires high computing power and is very time-consuming. Therefore, precise though it is, the Direct Numerical Simulation (DNS) method is rarely used.

The most economical are Reynolds Averaged Navier-Stokes (RANS) methods, which average vector and scalar quantities over time, introducing time mean values and value deviations [20], described by Eq. (4):

$v_i=\bar{v}_i+v_i^{\prime}$          (4)

They require the introduction of additional equations describing the turbulent stress tensor:

$\frac{\partial \rho}{\partial t}+\frac{\partial \rho}{\partial x_i}\left(\rho \bar{v}_l\right)=0$           (5)

$\left(\rho \bar{v}_l\right)+\frac{\partial \rho}{\partial x_j}\left(\rho \bar{v}_l \bar{v}_J\right)$ $=-\frac{\partial \bar{p}}{\partial x_i}+\frac{\partial}{\partial x_j}\left[\mu\left(\frac{\partial \bar{v}_l}{\partial x_j}+\frac{\partial \bar{v}_J}{\partial x_i}-\frac{2}{3} \delta_{i j} \frac{\partial \overline{v_k}}{\partial x_k}\right)\right]$$+\frac{\partial}{\partial x_j}\left(-\rho \overline{v_l^{\prime} v_J^{\prime}}\right)+\rho g_i+S_i$       (6)

Models based on the Boussinesq hypothesis [21] introduce the concept of turbulent viscosity to describe the turbulent stress tensor, defined as follows:

$-\rho \overline{\bar{v}_l^{\prime} v_J^{\prime}}=\mu_t\left(\frac{\partial \bar{v}_l}{\partial x_j}+\frac{\partial \bar{v}_j}{\partial x_i}\right)-\frac{2}{3}\left(\rho k+\mu_t \frac{\partial v_k}{\partial x_k}\right) \delta_{i j}$    (7)

The two most important Boussinesq models are k-$\varepsilon$ and k-w. In these models two other transport equations are added: one for the specific kinetic energy of turbulence k (which is the sum of the average kinetic energy and the instantaneous kinetic energy value, describing the turbulence fluctuations), and one for either the magnitude of the kinetic energy dissipation $\varepsilon$, or the specific dissipation rate w. The turbulent viscosity is computed as a function of k and $\varepsilon$, or k and w.

The standard k-$\varepsilon$ model [22] is regarded as relatively easy to use and it leads to stable computations but its disadvantages are the simplified description of energy dissipation and the fact that it only works for flows with developed turbulence, therefore its applicability for specific cases always needs to be verified. The k-w model, developed by Wilcox [20] possesses similar properties to k-$\varepsilon$, with better accuracy in lower Reynolds numbers. It is, however, highly sensitive to changes in the boundary conditions at the inlet with free flow.

This problem is solved by its modification, SST (Shear Stress Transport) k-w, proposed by Menter [23]. It is a combination of the k-$\varepsilon$ model, which is used in the fluid core, while k-w is applied in the boundary layer. Governing equations and standard parameters for this model are as follow:

$\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_i}\left(\rho k \bar{v}_i\right)$$=\frac{\partial}{\partial x_j}\left(\Gamma_k \frac{\partial k}{\partial x_j}\right)+G_k-Y_k+S_k+G_b$          (8)

$\frac{\partial}{\partial t}(\rho \omega)+\frac{\partial}{\partial x_i}\left(\rho \omega \bar{v}_i\right)$$=\frac{\partial}{\partial x_j}\left(\Gamma_\omega \frac{\partial \omega}{\partial x_j}\right)+G_\omega-Y_\omega+D_\omega+S_\omega+G_{\omega b}$              (9)

G_k and G_w are terms responsible for the production of k and w. $\Gamma_k$ and $\Gamma_\omega$ represent the effective diffusivity. Y_k and Y_w describe the dissipation of k and w due to turbulence. D_w is a cross-diffusion term. S_k and S_w are user-defined source terms. G_b and G_bw account for buoyancy.

$\Gamma_k=\mu+\frac{\mu_t}{\sigma_k}$    (10)

$\Gamma_\omega=\mu+\frac{\mu_t}{\sigma_\omega}$       (11)

$\sigma_k=\frac{1}{\frac{F_1}{\sigma_{k, 1}}+\frac{1-F_1}{\sigma_{k, 2}}}$        (12)

$\sigma_\omega=\frac{1}{\frac{F_1}{\sigma_{\omega, 1}}+\frac{1-F_1}{\sigma_{\omega, 2}}}$      (13)

where, $\quad \sigma_{k, 1}=1.176 ; \quad \sigma_{\omega, 1}=2.0 ; \quad \sigma_{k, 2}=1.0 ; \quad \sigma_{\omega, 2}=1.168 ; \quad F_1=$ $\tanh \left(\Phi_1^4\right) ; \quad \Phi_1=\min \left[\max \left(2 \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]$; $D_\omega^{+}=\max \left[2 \rho \frac{1}{\sigma_{\omega, 2}} \frac{1}{\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j}, 10^{-10}\right] ; D_\omega^{+}$is the positive portion of the cross-diffuse term; $y$ is the distance to the next surface.

$\mu_t=\frac{\rho k}{\omega} \frac{1}{\max \left[\frac{1}{\alpha^*}, \frac{S F_2}{a_1 \omega}\right]}$      (14)

where, $F_2=\tanh \left(\Phi_2^2\right) \Phi_2=\max \left[2 \frac{\sqrt{k}}{0.09 \omega y}, \frac{500 \mu}{\rho y^2 \omega}\right]$; a₁=0.31.

$\alpha^*=\alpha_{\infty}^*\left(\frac{\alpha_0^*+\frac{R e_t}{R_k}}{1+\frac{R e_t}{R_k}}\right)$         (15)

where, $R e_t=\frac{\rho k}{\mu \omega} ; R_k=6 ; \alpha_0^*=\frac{\beta_i}{3} ; \beta_i=0.072$.

For high-Reynolds number $\alpha^*=\alpha_{\infty}^*=1$.

$G_k=-\rho \overline{v_i^{\prime} v_j^{\prime}} \frac{\partial \bar{v}_j}{\partial x_i}=\mu_t S^2$          (16)

$S \equiv \sqrt{2 S_{i j} S_{i j}}$          (17)

$G_\omega=\frac{\alpha \alpha^*}{v_t} G_k$        (18)

$\alpha=\alpha_{\infty}\left(\frac{\alpha_0+\frac{R e_t}{R_\omega}}{1+\frac{R e_t}{R_\omega}}\right)$      (19)

where, R_ω=2.95; $\alpha_0=\frac{1}{9}$.

$\alpha_{\infty}=F_1 \alpha_{\infty, 1}+\left(1-F_1\right) \alpha_{\infty, 2}$          (20)

where, $\alpha_{\infty, 1}=\frac{\beta_{i, 1}}{\beta_{\infty}^*}-\frac{\kappa^2}{\sigma_{\omega, 1} \sqrt{\beta_{\infty}^*}} ; \alpha_{\infty, 2}=\frac{\beta_{i, 2}}{\beta_{\infty}^*}-\frac{\kappa^2}{\sigma_{\omega, 2} \sqrt{\beta_{\infty}^*}} ;$ k=0.41; β_i,1=0.075; β_i,2=0.0828.

$Y_k=\rho \beta^* k \omega$       (21)

where, $\beta^*=\beta_i^*\left[1+\zeta^* F\left(M_t\right)\right] ; \beta_i^*=\beta_{\infty}^*\left(\frac{\frac{4}{15}+\left(\frac{R e_t}{R_\beta}\right)^4}{1+\left(\frac{R e_t}{R_\beta}\right)^4}\right) ; \zeta^*=$ $1.5 ; R_\beta=8 ; \beta_{\infty}^*=0.09$.

F(M_t) is the compressibility function. In the incompressible flows $\beta^*=\beta_i^*$. In the high-Reynolds number form $\beta_i^*=\beta_{\infty}^*$.

$Y_\omega=\rho \beta \omega^2$         (22)

$\beta=\beta_i\left[1-\frac{\beta_i^*}{\beta_i} \zeta^* F\left(M_t\right)\right]$           (23)

$\beta_i=F_1 \beta_{i, 1}+\left(1-F_1\right) \beta_{i, 2}$         (24)

$G_{\omega b}=\frac{\omega}{k}\left[(1+\alpha) C_{3 \varepsilon} G_b-G_b\right]$          (25)

The SST k-w model is able to accurately represent both laminar and turbulent flows [24].

Menter et al. [25] developed the Intermittency Transition Model, specifically to describe the laminar-to-turbulent transition. The model is later coupled with SST k-w. Added governing equations are as follow:

$\frac{\partial(\rho \gamma)}{\partial t}+\frac{\partial\left(\rho U_j \gamma\right)}{\partial x_j}=P_\gamma-E_\gamma+\frac{\partial}{\partial x_j}\left[\left(\mu+\frac{\mu_t}{\sigma_\gamma}\right) \frac{\partial \gamma}{\partial x_j}\right]$          (26)

$P_\gamma=F_{\text {length }} \rho S \gamma(1-\gamma) F_{\text {onset }}$        (27)

where, S is a strain rate magnitude and F_length=100.

$E_\gamma=c_{a 2} \rho \Omega_\gamma F_{t u r b}\left(c_{e 2} \gamma-1\right)$         (28)

where, $c_{a 2}=0.06 ; c_{e 2}=50 ; \Omega_\gamma=1.0 ; F_{\text {onset }}=\max \left(F_{\text {onset2 }}-F_{\text {onset3 }}, 0\right)$; $F_{t u r b}=e^{-\left(\frac{R_T}{2}\right)^4} ; R_T=\frac{\rho k}{\mu \omega} ; R e_V=\frac{\rho d_w^2 S}{\mu}$.

To calculate the critical momentum thickness Reynolds number, the freestream turbulence intensity and pressure gradient are approximated locally:

$R e_{\theta c}=f\left(T u_L, \lambda_{\theta L}\right)$

$R e_{\theta c}\left(T u_L, \lambda_{\theta L}\right)=C_{T U 1}$

$+C_{T U 2} \exp \left[-C_{T U 3} T u_L F_{P G}\left(\lambda_{\theta L}\right)\right]$

$T u_L=\min \left(100 \frac{\sqrt{\frac{2 k}{3}}}{\omega d_w}, 100\right)$         (29)

$\lambda_{\theta L}=-7.57 \cdot 10^{-3} \frac{d V}{d y} \frac{d_w^2}{v}+0.0128$

To ensure numerical robustness, $\lambda_{\theta L}$ is bounded:

$\lambda_{\theta L}=\min \left(\max \left(\lambda_{\theta L},-1.0\right), 1.0\right)$        (30)

$F_{P G}\left(\lambda_{\theta L}\right)$$=\left\{\begin{array}{r}\min \left(1+C_{P G 1} \lambda_{\theta L}, C_{P G 1}^{l i m}\right), \lambda_{\theta L} \geq 0 \\ \min \left(1+C_{P G 2} \lambda_{\theta L}+C_{P G 3} \min \left[\lambda_{\theta L}+0.0681,0\right], C_{P G 2}^{l i m}\right), \lambda_{\theta L}<0\end{array}\right.$           (31)

$F_{P G}=\max \left(F_{P G}, 0\right)$           (32)

Model constants: $\quad C_{T U 1}=100.0 ; \quad C_{T U 2}=1000.0 ; \quad C_{T U 3}=1.0 ;$ $C_{P G 1}=14.68 ; C_{P G 2}=-7.34 ; C_{P G 3}=0.0 ; C_{P G 1}^{l i m}=1.5 ; C_{P G 2}^{l i m}=3.0$.

The model was established for the external flows, to adjust it for internal flows Abraham et al. [26] proposed the change to constant values as follows: C_TU1=75.0; C_TU3=1.37; β*=0.092.

All described models are suitable for modelling transitional flow, that occurred during the experimental research (see Section 3).

2.2 Porous media model

The porous media model is based on the Darcy-Forchheimer equation, which describes fluid flow through porous media that differs from Darcy's law. At a high Reynolds number, the inertial forces grow more prominent, and the relationship between pressure gradient and seepage velocity becomes non-linear. When those forces are of higher order than viscous forces, non-Darcy flow can be described with the empirical Forchheimer equation [27]:

$-\nabla P=\frac{\mu}{\alpha_p} v+\beta_p \rho v^2$     (33)

The inertial $\left(\beta_p\right)$ and viscous (α_p^-1) coefficients of the porous media model need to be calculated from experimental data.

The porous media model in ANSYS Fluent adds the momentum sink to the governing momentum equations. The momentum source term consists of a viscous loss (Darcy) term and an inertial loss term, as presented in the following equation [17]:

$S_i=-\left(\sum_{j=1}^3 D_{i j} \mu v_j+\sum_{j=1}^3 C_{i j} \frac{1}{2} \rho|v| v_j\right)$         (34)

By default, the built-in model uses superficial velocity, so the porous media coefficients from the experimental data need to be calculated accordingly. The effect of the porous zone on the turbulence field is only approximated. In the default approach, the solid medium has no effect on turbulence generation, which is a reasonable assumption if its permeability is large, i.e., typical of turbulators.

The CFD simulations were performed using the ANSYS software package, particularly the Fluent solver. Only the inner tube of the heat exchanger was modelled since this is where the turbulators were located and the pressure was measured.

4.1 Selection of a numerical mesh

The first stage of the simulation is to create a numerical mesh. The mesh independence study ensures that mesh quality does not affect results. Examples of generated meshes and their properties are shown in Table 2. Polyhedral meshes allow good geometry fit and high quality while significantly reducing the number of cells. Since a high gradient of velocity values was expected on the walls, each time the boundary layer of ten cells was added. To test the mesh two turbulators were chosen for the calculations – one representing turbulators affecting mainly the boundary layer, and one taking up the whole volume of the pipe. Simulations were performed for each mesh by the SST k-w turbulence model. The relative errors of the static pressure drop were calculated for each data point.

$\delta_p=\frac{\Delta p_{c a l c}-\Delta p_{e m p}}{\Delta p_{e m p}} \cdot 100 \%$     (35)

The results are shown in Figure 4. The medium-size mesh, selected due to the limited computational power available, was chosen to find out the balance between the calculation time and accuracy of the results. The average relative error drops only by less than 1 percentage point for the larger mesh, while the calculations lasted over 4 times longer.

Table 2. Properties of the generated meshes

Numerical mesh	2.1.png	2.2.png	2.3.png
Number of cells	127 thousand	1.52 million	5.38 million
Maximum size of a surface cell	1.46 mm	0.46 mm	0.26 mm
Max Skewness	0.36	0.47	0.55
Average Skewness	0.02	0.01	0.01
Average relative error for the brush turbulator with Pall Rings	10.7%	9.1%	9.0%
Average Relative Error for the spring turbulator – 8 mm diameter, 2 mm pitch	21.9%	19.4%	18.6%

4.png

Figure 4. Static pressure drop for spring turbulator ϕ8 mm (ST) and brush turbulator with Pall rings (BT) for various mesh density

4.2 Selection of a turbulence model

With the mesh being selected, three relevant turbulence models offered by ANSYS Fluent, suitable for simulating transitional flows (see Section 2.1), were compared. The average relative errors along with their standard deviations for each model are presented in Table 3. The standard deviation indicates how close individual errors are to the reported average.

It can be observed that altering the constants in the intermittency model does not affect the results at all in the discussed flow regime. Furthermore, in our case, the standard SST k-w model performs just as well, while reducing computational effort (no need to solve additional equations). This model is recommended by ANSYS. It was chosen as the most economical option for future calculations, given that the calculation time was shorter, and the difference in the results negligible.

Table 3. Average relative error and standard deviation for each turbulence model for the numerical mesh of 1.52 million cells for different turbulators

Turbulence models: 1 - SST k-w; 2 - Intermittency model; 3 - Intermittency model with changed constants

4.3 Simulations of the pressure drop for turbulators

The coefficients presented in Table 1 were implemented in the porous medium model and simulations were performed. The results of the average relative error and standard deviation for each turbulator are presented in Table 4.

Table 4. Average relative error and standard deviation for each turbulator for the SST k-omega turbulence model

The degree to which the porous medium model approximates the pressure drop in the turbulators can be correlated with the values of the inertial coefficients presented in Table 1. The greater the influence of inertial forces, the more the turbulator resembles the porous media, and better accuracy is achieved. Twisted tape turbulators have the lowest inertial coefficients. They do not block the fluid flow in the way a porous media would; hence the low accuracy. Stamped sheet and spring turbulators have similar relative errors and medium inertial coefficients. Pall rings could resemble porous media of medium porosity and the relative error corresponds to that. Brush turbulators occupy most of the space, have the highest inertial coefficients, and the lowest errors. It can be concluded that turbulators that mainly affect turbulence by creating an obstacle for the flow at the fluid core are better described by the porous media model than turbulators that act on the boundary layer.

This research provided a baseline understanding of the problem of simulating turbulators using the porous media model, which has not been discussed so far. It is a good starting point for future exploration.

The porous media model works best with turbulators described by high inertial coefficients. Their common trait is that they all occupy a high percentage of the cross-section of the pipe and cause the highest pressure drop.

The porous media model cannot be used to describe twisted tape turbulators. Their effectiveness lies mainly in directing the flow towards the walls of the device and breaking the boundary layer, which is not modelled properly.

Modelling the real geometry of the turbulators is very time-consuming, and the results are highly dependent on the quality of the mesh; therefore, it would be advisable to simplify the process.

It is worth noting that the porous media model returns a higher value of the pressure drop than in reality, which is fine for engineering applications, since continuous flow must always be ensured while designing an installation, and it is better to overestimate than to underestimate the pressure drop.

Modelling the pressure drop in turbulators using the porous media approach gives promising results and should be further developed. The next stage of the research would be to introduce heat exchange into the simulation.