Numerical method

A fully developed channel flow is simulated in the domain shown in Figure 1, where the boundaries on X and Z direction are periodic boundaries and the boundary on Y direction is no-slip isothermal wall boundary - mean flow is in the X direction. The scale of the domain is Lx×Ly×Lz=2πh×2h×πh, where h is the characteristic length of the channel.

Figure 1.

Calculation domain of the channel flow simulation.

This work is based on the solution of the three-dimensional Navier-Stokes (NS) equations (Equations 2–4), which are closed by the equation of state p=p(ρ,T). These equations are non-dimensionalised by half-channel height, h, globally averaged density, ρG≡(1/h)∫0hρ¯dy, and global averaged velocity, uG≡(1/hρG)∫0hρu¯dy, before being calculated.

where:

Non-dimensional parameters can be declared as:

For a real gas, all the thermodynamic parameters are two-variable functions of both temperature and density. The global density and wall temperature (ρG, Tw) are taken to determinate the global dynamic viscosity, μG, and heat capacity Cp,G during the non-dimensionalisation. The global Reynolds number is defined by ReG≡(ρGuGh/μG) and the global Prandtl number is defined by PrG≡(Cp,GμG/λG). Subsequently, expression of the non-dimensional NS equations is the same as the original expression.

For the fully developed channel flow there is a periodic boundary instead of inlet and outlet boundaries, so the pressure drop along the stream is balanced by an artificial added body force fi, which is equivalent to the mean pressure gradient in physical flow (Huang et al., 1995). The body force fi is non-zero only for the mean flow direction (i=1), and it is uniform in the domain. The body force is calculated in each time step by keeping the global averaged mass flow rate invariant in the calculation.

The DNS solver is developed based on a compressible DNS code from X. Li (Li et al., 2001) with an addition of a real-gas model. All the thermal properties in this model are two-variable functions of the local density and temperature, and they are calculated by the state-of-art Helmholtz energy EoS in REFPROP (Lemmon et al., 2007). Organic fluid cases are calculated by this real-gas model while the air case is calculated by a perfect-gas model. The NS equations are solved after non-dimensionalisation. For spectral discretisation, the convective flux derivatives are calculated by seventh-order upwind finite-difference scheme, and the viscous flux derivatives are calculated by sixth-order central scheme. For time advancement, it is solved by a third-order Runge-Kutta method. The validation of the simulation method is shown in the appendix.

The number of nodes in this work is Nx×Ny×Nz=128×257×128, and the overall number of nodes is 4.2 × 10⁶. The spatial resolution is chosen as DNS requirements. Zonta (Zonta et al., 2012) advised that (Δx/η)max≈12, (Δy/η)max≈2 and (Δz/η)max≈6, while Lee (Lee et al., 2013) advised that (Δx/η)max≈12.4, (Δy/η)max≈3.0 and (Δz/η)max≈7.9, where η is the smallest length scale in the flow. η is equal to the Kolmogorov length scale ηk when Prandtl number is lower than 1, while it equal to ηθ≈ηk/Pr when Prandtl number is higher than 1 (Monin et al., 1975). The spatial resolution of cases in this work is shown in Table 2, and it can fulfill the DNS requirement. The y+ of the first layer from wall boundary is below 0.5, and the first 10 layers are within y+=5.

Figure 2a and 2b show two-point correlation of velocity u and density ρ between points in distance Δx∈[0,Lx/2] at y+≈10, and Figure 2c and 2d show the two-point correlation between points in distance Δz∈[0,Lz/2]. It shows that all the two-point correlations are close to zero at Lx/2 and Lz/2, which means that the domain size Lx and Lz are large enough for the channel flow cases with periodic boundary.

Figure 2.

Two-point correlation of velocity u and density in streamwise and spanwise. correlation in streamwise (a) and (b); correlation in spanwise (c) and (d).

This work is solved by the High Performance Computer in Imperial College London. Each case uses 512 CPUs, and it takes 10–20 days to get converged. After that, the statistic is taken every Δt^=2 within a time range Δt^=1000.

Case set-up

There are no standard inlet and outlet boundaries in this domain for the periodic boundary condition, so the global averaged density ρG and mass flow rate m˙G=uGhρG are set as an alternative. The global averaged velocity is given as uG=MaGcG, where MaG is the global Mach number. cG is a reference global sound speed defined by wall temperature and global averaged density cG=crg(ρG,Tw) for organic fluids, and only wall temperature cG=cpg(Tw) for air. The dynamic viscosity μG is defined in the same fashion - μG=μrg(ρG,Tw) for organic fluids and μG=μpg(Tw) for air. μrg is calculated by the real-gas model, and μpg is calculated by the Sutherland viscosity law (Sutherland, 1893) shown in Equation 16, where T_s = 110.56 K is the effective temperature.

The initial condition is a turbulent flow field, in which the global averaged density ρG and mass flow rate m˙G are set. During the calculating, ρG and m˙G will not change. After each time step, ρG is calculated and compared with its value in the previous step. If there is a gap between them (caused by numerical error), the density at all points will plus or minus a same small density source term Δρ so as to keep ρG invariant after this step. m˙G is also calculated in each step and compared with its value in the previous step. If they are different, the body force will be changed so as to maintain an invariant m˙G.

The global Mach number MaG≡(uG/cG) and global Reynolds number ReG≡(ρGuGh/μG) are set at 2.25 and 4,880 for all seven cases. The air case follows the perfect-gas Equation of State, requiring only one thermodynamic condition - the wall temperature, T_w = 288 K. The two thermodynamic conditions chosen for the real-gas cases were wall temperature and global density. From R1 to R3, the set-up point moves closer to the supercritical region, shown by the temperature-entropy diagram in Figure 3. The red and blue lines in Figure 3 are the averaged thermodynamic property distribution in each real-gas case. Cv/R for each case is list in Table 3, and the value for MDM is the largest whilst air is the smallest. As the set-up point closer to the supercritical region Cv/R increases, although this change is small in comparison to the difference in values between R1233zd(E) and MDM. The fundamental derivative Γ shows that the minimum value is at R2 and M2 for the two organic fluids, respectively.

Figure 3.

Temperature-entropy diagram of R1233zd(E) and MDM. Dots refer to set-up cases, whilst lines refer to mean thermal property distribution. (a) R1233zd(E) (b) MDM.

Averaged velocity and thermodynamic properties

The variables are averaged on each y level, since the x and z directions are periodic. For analysis, the wall shear stress is defined as the averaged value τw≡[μ(∂u/∂y)]¯w and wall friction velocity is defined as uτ≡τw/ρ¯w. Huang et al (Huang et al., 1995) proposed a semi-local friction velocity uτ∗≡τw/ρ¯ as well, which takes compressibility effects into consideration.

Some important global results of each case are listed in Table 4. Bulk Reynolds number and bulk Mach number are defined as ReB≡ρGuGh/μ¯w and MaB≡uG/c¯w (Sciacovelli et al., 2017), and they are different from ReG and MaG in real gases (μ¯w≠μG).

The centre-line Reynolds number and Mach number are defined as ReC≡ρ¯Cu¯Ch/μ¯C and MaC≡uC/cC¯.The values show that ReB decreases as the set-up point is closer to the supercritical region (from R1 to R3 and from M1 to M3), while ReC increases. Both organic fluids have smaller ReB and larger ReC relative to the air baseline. The Mach number MaB firstly increases and then decreases as the set-up point is closer to the supercritical region (from R1 to R3 and from M1 to M3), while the centreline Mach number shows minimal change across cases.

The friction Reynolds number is defined as Reτ≡ρ¯wuτh/μ¯w. The friction Reynolds number is lower in the organic fluid cases than air, while lowest in the supercritical cases (R3 and M3). The semi-local friction Reynold number Reτ∗≡ρ¯uτ∗h/μ¯ has a higher value for the two organic gases than air, and Reτ∗ is highest in the supercritical cases R3 and M3 for each fluid.

Reynolds averaged velocity and distribution versus non-dimensional wall distance y+≡ρ¯wuτy/μ¯w is shown in Figure 4a The position y+=0 is the wall, while the maximum y+ is the central layer of the channel. It shows that the organic fluid cases are closer to the log law line u+=0.25 ln (y+)+5.5 in the outer layer than air. In each fluid, as the case moves closer to the supercritical region, the velocity distribution increases in deviation from the log law line. After the application of a Van Driest transformation (Equation 17) (Van Driest, 1951) good agreement to the log law in the outer layer is observed (shown in Figure 4b).

Figure 4.

Distribution of x-direction averaged velocity along wall distance. black dashed lines refer to u+=y+ in the viscous sublayer and u+=0.25ln⁡(y+)+5.5 in the outer layer.

The averaged temperature, density and pressure along the wall-normal direction are shown in Figure 5, where the abscissa is the non-dimensional wall distance y∗≡ρ¯uτ∗y/μ¯. Figure 5a shows that the temperature in the central region of the channel is higher than wall temperature, and centre-to-wall temperature ratio T¯C/Tw is lower in the organic fluid cases than in the air case (shown in Table 4).

Figure 5.

Distribution of averaged temperature, density and pressure along the wall-normal direction. (a) Averaged temperature. (b) Averaged temperature (zoomed). (c) Averaged density. (d) Averaged pressure.

For each organic fluid, T¯C/Tw decreases as the point moves closer to the supercritical region (R1 to R3 and M1 to M3). T¯C/Tw correlated positively to PrEc in these cases (Prandtl number Pr≡Cpμ/λ, Eckert number Ec≡c2/CpT), which means the larger PrEc the larger T¯c/Tw will be.

When averaging and integrating the energy equation, the heat transfer as a function of wall distance can be represented by Equation 18 [12]. In the expression, Fu(y) is a velocity distribution function which equal to 0 at y^=0 and equal to 1 at y^=1. Since the heat conductivity fluctuation is small compared with averaged conductivity, the gradient of averaged temperature can be derived to Equation 19.

By assuming the velocity distribution function and the thermal conductivity ratio the same in all these cases, the non-dimensional temperature gradient is approximately proportional to the global Prandtl number multiplied by global Eckert number PrGEcG whose values are shown in Table 4. For a perfect gas, EcG=γR/Cp=γ−1. The organic fluids have higher molecular complexity, leading to a larger Cp/R and therefore a smaller Ec. For each organic fluid, Pr and Ec are variable, and the closer to the super critical region the smaller the value of PrEc.

The density is higher at the channel wall than the channel centre (shown in Figure 5c, and the value of the wall-to-centre density ratio ρ¯w/ρ¯C is shown in Table 4. The largest density ratio also appears in the air case, which is the same for the temperature ratio mentioned above. However, for each organic fluid, the density ratio increases as the set-up point moves closer to the supercritical region, while the temperature ratio has a reverse trend.

Since the pressure is nearly constant (p¯/p¯w>0.985 shown in Figure 5d), the density gradient is approximately equal to the temperature gradient multiplied by −(∂ρ^/∂T^)p (Equation 21). The value of −(∂ρ^/∂T^)p is shown in Table 4. It is equal to 1 in perfect gas, but it increases significantly as the set-up point moves closer to the supercritical region in real gases(from R1 to R3, and M1 to M3). The effect of the increase in −(∂ρ^/∂T^)p is larger than that of the decrease in temperature gradient, so the density drop is larger in cases closer to the supercritical region.

The pressure ratio has a minimum value at y∗≈70 for all the cases (shown in Figure 5d). As the set-up point moves closer to the supercritical region, the maximum pressure drop firstly increases and then decreases in R1233zd(E) (from R1 to R3), while it increases monotonically in MDM (from M1 to M3). This trend of maximum pressure drop is nearly negatively related to the averaged pressure at wall p^¯w (shown in Table 3), and this relation can be explained by the Reynolds averaged momentum equation in the y direction (Equation 22) in the following paragraph.

After integrating Equation 22 along wall distance y, the expression of for pressure ratio shown in Equation 23 is found. The pressure ratio is associated to 1/p^¯w and ∫0y(d(ρ¯v′′v′′~)/dy)dy. If the later part differs slightly between cases, it can be concluded that the larger the wall pressure p^¯w is the smaller the pressure drop will be.

The distribution of averaged dynamic viscosity ratio μ¯/μ¯w is shown in Figure 6. It shows that the viscosity increases as y∗ increases in the Air baseline, R1 and M1, while it decreases in Case R3 and M3. Sciacovelli et al (Sciacovelli et al., 2017) proposed that the former trend is a gas-like behaviour while the latter is a liquid-like behaviour, and that liquid-like behaviours can be observed in dense-gases.

Figure 6.

Distribution of averaged dynamic viscosity along the wall-normal direction.

From the result of Figure 6, it shows that the liquid-like behaviour only exists in the supercritical cases (R3 and M3) of the two organic fluids. By implication there is a transition line from gas-like behaviour to liquid-like behaviour in each organic fluid, which is (∂μ/∂T)p=0 shown in Figure 7. The slope of the transition line is larger than that of the isobar lines. When fluid works below the transition line, the dynamic viscosity increases with wall distance; otherwise when fluid works above that line, the dynamic viscosity decreases as wall distance is increased.

Figure 7.

Transition line for R1233zd(E) and MDM. (a) R1233zd(E). (b) MDM.

Reynolds stress

The distributions of Reynolds stresses u′′u′′~, v′′v′′~, w′′w′′~ and u′′v′′~ normalized by uτ∗2 are shown in Figure 8. The normal Reynolds stress in the streamwise direction, u′′u′′~, shows a slight deviation (less than 4%) between cases (shown in Figure 8a), while the normal Reynolds stress in spanwise and wall-normal directions (v′′v′′~and w′′w′′~) show larger values of deviation (Figure 8b and 8c).

Figure 8.

Distribution of reynolds stresses.

The peak values of v′′v′′~ in real-gas cases are larger than that in air case by 7–20%. For each fluid, when the real-gas effect increases(from R1 to R3, or M1 to M3), the peak value will increase. Meanwhile, this impact is stronger in R1233zd(E) than in MDM. Real-gas effect also influence the peak values of w′′w′′~, which are larger by 10–21% in real-gas cases. The trend of this impact is the same as that of the real-gas effect on v′′v′′~.

The impact of real-gas effect on the normal Reynolds stress in spanwise and wall-normal directions (v′′v′′~ and w′′w′′~) is associated with the viscosity profile and bulk Reynolds number ReB. Patel et al. (Patel et al., 2015) proposed that if the centre-to-wall viscosity ratio less than 1, the normal Reynolds stress in spanwise and wall-normal directions will increase, and Sciacovelli et al (Sciacovelli et al., 2017) proposed that the increase of bulk Reynolds number ReB will lead to a larger peak value of normal Reynolds stress. Both these two effects influence the normal Reynolds stress in the cases of this work, while the viscosity ratio effect seems more predominant in this work, since supercritical case R3 (with the lowest μ¯c/μ¯w and lowest ReB) has the largest maximum v′′v′′~ and w′′w′′~ among all cases.

The Reynolds shear stress u′′v′′~ is larger in two organic fluids by (5–8%) in comparison to air, while the difference between organic fluid cases is not significant. The maximum value of Reynolds shear stress u′′v′′~ is associated with the viscous shear stress. According to the Reynolds averaged momentum equation in the y direction (Equation 24), the viscous shear stress and Reynolds shear stress balance the body force together. Since the dynamic viscosity for real-gas cases is smaller in the outer layer, which leads to a smaller viscous shear stress, so the Reynolds shear stress becomes larger compared with air case.

Turbulence energy

The Reynolds averaged turbulence energy equation for fully developed channel flow is shown in Equation 25. The first term is turbulence energy production Pk; the second term is the summary of turbulent diffusion, velocity-pressure diffusion, and viscous diffusion; the third term is energy dissipation ϵk; and the last three terms are compressibility-related terms. The turbulence energy production term Pk and dissipation term ϵk are discussed in this work, which are defined in Equations 26 and 27.

The distribution of turbulence energy production Pk is shown in Figure 9a. It shows that the fluid property do not significantly influence non-dimensional turbulence energy production, with a maximum Pk difference of less than 5% across all cases. The distribution of turbulence energy dissipation ϵk is shown in Figure 9b. The difference between the air and real-gas cases are significant in the viscous sublayer and buffer layer (y∗<30), with the magnitude decreasing in the outer layer (y∗>30). The non-dimensional dissipation at the wall is 5–13% higher in real-gas cases than air. For each organic case, the non-dimensional dissipation in the viscous sublayer and buffer layer increases as the real-gas effect increases (from R1 to R3; M1 to M3).

Figure 9.

Distribution of turbulence energy (a) Production, (b) Dissipation.

The turbulence energy production-to-dissipation ratio −Pk/ϵ−1 is shown in Figure 10. A clear peak is demonstrated at y∗≈10 across all cases, followed by a reduction to nearly zero as y∗≈50 is reached. In the region y∗>50, the air case remains stable before reducing to −1, while those for real-gas cases will increase slowly before demonstrating the same decrease (albeit at higher values of y∗). The position of the start point for the second drop is related to Reτ,C∗, and the larger the Reτ,C∗ the later the second drop will be.

Figure 10.

Distribution of turbulence energy production - to - dissipation ratio.

Turbulence viscosity hypothesis

The evaluation of the turbulence viscosity hypothesis for the standard two-function RANS model (k−ϵ and k−ω) is a focus of study in this work. The turbulence viscosity νt is defined as the ratio of Reynolds shear stress to the gradient of the averaged velocity, whose simplified expression in the full developed channel flow is shown in Equation 28. In the standard k−ϵ model (Launder and Spalding, n.d.), the turbulence viscosity νt is calculated by the turbulence kinetic energy k and dissipation ϵ (Equation 29), which are solved by equations for k and ϵ. The constant Cμ is equal to 0.09.

The ratio νtϵ/k2 is calculated for each case and shown in Figure 11. It shows that in the region y∗<70, the ratio is lower than 0.09 for all cases, but approximately equal to 0.09 in the region y∗>70 within deviation bands of ±10%.

Figure 11.

Distribution of νtϵk/kt2 across all cases.

With the data for R1233zd(E) and MDM collected so tightly around the constant Cμ line across all cases, there is little evidence to suggest that real-gas effects have a significant impact on this value. A conclusion can therefore be made that the assumption Cμ=0.09 is still suitable in the standard two-function RANS model (k−ϵ and k−ω) within the outer layer region at y∗>70 for the organic fluids within this study.

Greeks

Subscripts

Accents