CFD Applications in Ground Source Heat Pump System

Abstract

In ground source heat pump (GSHP) system, computational fluid dynamics (CFD) is commonly used to conduct simulation analysis of its operating characteristics. Particularly, ground heat exchanger (GHE) is the most core component of GSHP system, and the heat transfer characteristics of which with soil around will directly affect the efficiency of the entire system. Thus, CFD is always applied to predict the process of heat transfer around GHE and its influence on heat exchange process. In this chapter, a 3-D numerical model considering dynamic surface condition and initial soil temperature distribution is developed to investigate the thermal performance of helix ground heat exchanger (HGHE) on basis of CFD, and the main influencing factor (inlet water temperature) is studied with the established model. In addition, the experimental investigation is carried out to verify the accuracy of the model. The results are of great significance for exploring the application of CFD in GSHP system.

1. Introduction

Shallow geothermal energy resource is a kind of renewable resources, mainly referring to the low-temperature heat energy, also known as geothermal energy, in the earth’s shallow surface within hundreds of meters (less than 200 m). Different from the traditional deep geothermal energy, it refers to the stratum thermal energy in the general layer, the temperature of which is less than 25°C, and its energy mainly comes from solar radiation and earth gradient warming. Compared with deep geothermal energy in the traditional sense, shallow geothermal energy is not affected by geological factors and generally exists in the surface of the earth, with many advantages such as wide distribution and renewable [1, 2, 3].

Due to the low temperature of shallow geothermal energy, which is not easy to be extracted, it was not used by people in the early years. With the progress of science and technology and the attention paid to the impact of natural environment, as a renewable, clean, and huge kind of energy, it has gradually been widely paid attention to. At present, the shallow geothermal energy is mainly used in the air conditioning system of buildings. The method is to extract the low-grade underground shallow geothermal energy through the GSHP system and make use of it. Figure 1 shows the schematic diagram of the GSHP system. By burying the GHE underground, it makes heat exchange with the relatively constant temperature stratum. In winter, it takes the earth as the high-temperature potential heat source and uses a small amount of high-grade energy (electric energy) to extract heat from the earth and provides heat for the user. On the contrary, in summer, it takes the earth as the low-temperature potential heat source and releases the excess heat of the user into the earth [4, 5, 6, 7, 8].

Figure 1.

Schematic of ground source heat pump system.

Traditional GSHP system can be divided into two kinds: water source heat pump technology and soil source heat pump technology. The heat source of the water source heat pump system is the underground or surface water, and the vast majority of the water source heat pump system directly extracts underground or surface water, which is an open system. The soil source heat pump system is a closed system by burying the GHE in the ground and exchanging heat with the surrounding rock and soil through the flow of the working medium inside the GHE. Because the underground or surface water used by water source heat pump system contains microorganisms, calcium and magnesium ions, and sulfate ions, the scale and blockage always occur in the system, resulting in poor system stability. At the same time, it is also restricted by hydrogeological conditions; thus, for the sake of realizing the development and utilization of shallow geothermal energy resources, the soil source heat pump system is widely used in the project at present.

Buried GHE is the most core component in the GSHP system. As shown in Figure 2a, the application research of the traditional U-type GHE started earlier and developed more mature. Therefore, the U type GHE is the most widely used, and most GSHP projects adopt this type of buried pipe. Compared with the U type GHE, the HGHE is a new type of geothermal exchanger, as shown in Figure 2b. Because the spiral pipe is arranged along the cylindrical wall, its heat transfer area and heat transfer amount are large, with low initial investment, high heat transfer efficiency, and other significant advantages, and HGHE with its unique advantages has been rapidly developed in recent years, effectively solving the economic and technical problems of traditional U-type GHE.

Figure 2.

Schematic of ground source heat pump system. (a) Traditional U-type GHE, (b) New-type HGHE.

Due to the complicated hydrogeological conditions involved in the actual engineering of GHE, it is difficult to obtain the heat transfer process of buried pipe by experimental means due to the limitation of the region and testing conditions. For this reason, scholars around the world have carried out a lot of numerical heat transfer simulation works on GHE, mainly by using CFD technology to achieve the simulation research on the heat transfer characteristics and heat transfer process of buried pipes, for the sake of maximizing the heat transfer efficiency and improving the overall efficiency of the GSHP system. Angelo Zarrella [9, 10, 11] et al. established a numerical heat transfer model based on U-type and helix GHEs by using the electrical resistance analysis method. They considered the influence of external environment on soil boundary conditions and conducted experimental research. The model mainly considers the convective heat transfer of helical fluid, and the radial and axial heat conduction of soil and pile foundation, and ignores the heat transfer along the Angle direction. Salsuwanda Selamat [12] aiming at the design optimization problem of horizontal HGHE, established a three-dimensional numerical model and used CFD software for simulation calculation to simulate and optimize each parameter. Gyu-Hyun Go et al. [13] established a three-dimensional numerical model of horizontal HGHE and simulated its heat transfer characteristics. The accuracy of the model was verified through indoor thermal response experiment. Qiang Zhao et al. [14] established a three-dimensional unsteady heat transfer model reflecting the heat transfer process of fluid and soil in the tube for the GHE and used the finite element method to solve and calculate. Then, the heat transfer characteristics of U-type, W-type, and helix GHEs are simulated and calculated, respectively. The research shows that the heat transfer efficiency of HGHE is higher than that of other two types of heat exchangers.

In order to further elaborate and introduce the application of CFD technology in the GSHP system, this chapter takes the new-type HGHE as an example and establishes its three-dimensional numerical model, which fully considers the external dynamic environmental conditions and initial soil temperature conditions. Then, the ANSYS_FLUENT software, which is commonly used in CFD technology, is used to simulate and study the heat transfer characteristics of HGHE under different inlet water temperature conditions. Finally, the accuracy of the numerical model is verified by experimental means.

2. Description of numerical model

2.1 Physical model

The research object of this chapter is made of polyethylene (PE) pipe, which is applied in Chongqing, China. The construction process is strictly in accordance with the relevant requirements of Technical Specifications for Ground Source Heat Pump System Engineering. And the specific parameters are shown in Table 1.

Parameters	Spiral pitch in the depth direction/m	Diameter/m	Height/m	Depth of the top of HGHE from ground surface/m	Inside diameter of coil pipe/m
HGHE	0.08	1.14	2.00	1.00	0.02

Table 1.

Main structural parameters of HGHE.

Using the Design Modeler module in ANSYS_FLUENT, which is a professional finite element analysis tool, this work established the geometric model of HGHE. According to the real size, the surrounding soil was set as a cylinder in the physical model. And the geometric model is shown in Figure 3a.

Figure 3.

Mesh division of physical model of HGHE. (a) Physical model of HGHE. (b) Mesh division of HGHE.

Then, the meshing module of ANSYS_FLUENT software is used to divide the grid of the built geometric model of HGHE. The grid division follows the following principles: If the temperature field or velocity field in a certain position and direction has significant fluctuations, the area should be divided into dense grids. On the contrary, if the fluctuation of temperature field or velocity field is not obvious, the meshing of this area can be loosened appropriately. Because the scale of fluid flow and heat transfer in the spiral buried pipe are small and fluctuate greatly during heat transfer, the grid of the spiral buried pipe and its nearby soil area are encrypted, as shown in Figure 3b. The specific grid division is as follows:

1. Hexahedral mesh is selected for mesh division of HGHE. Hexahedral mesh can reduce the number of elements and improve the convergence speed and analysis accuracy. Because the flow direction of the circulating medium in the tube changes at all times, the flow field fluctuates significantly, and local encryption is carried out during grid division to avoid large inclination angle of the grid in this area, which will promote the convergence speed and the accuracy of the calculation results.

2. The soil around the HGHE is hexahedral mesh, in which the area adjacent to the buried pipe is locally encrypted to avoid large inclination angle of the grid in this area, which will also promote the convergence speed and the accuracy of the calculation results.

3. Mesh with higher quality can accelerate the convergence speed, reduce the error in solving, avoid the divergence of numerical simulation results, and ensure the accuracy of results. ANSYS_FLUENT software is required to ensure that the minimum orthogonal quality is greater than 0.1 or the maximum inclination is less than 0.95. Most of the grids in the model established in this chapter have good quality, and a few grids with poor quality are in the acceptable range of ANSYS_FLUENT. Therefore, the grids used for numerical calculation can meet the requirements of ANSYS_FLUENT.

2.2 Model assumption

To simplify the solution, the following assumptions are made for the numerical model:

1. The soil around the buried pipe is regarded as an idealized uniform porous medium, and the thermal physical parameters of the soil remain constant in the heat transfer process.

2. The solid skeleton of the soil has no deformation, and the inner space is interconnected, with uniform porosity.

3. The contact thermal resistance between buried pipe and soil is not considered.

4. The influence of soil seepage and water migration is not considered.

5. Ignoring the influence of the HGHE around, the adiabatic boundary conditions were assumed around the cylindrical soil.

6. The local transient thermal equilibrium state is considered between soil and buried pipe.

2.3 Governing equations

The soil around the HGHE is regarded as a uniform porous medium, while the flow state of the medium in the spiral buried pipe is turbulent. Therefore, the governing equations of fluid flow and heat transfer in pipe and heat transfer of soil are established respectively.

2.3.1 Governing equations of fluid flow and heat transfer in pipe

Due to the turbulent flow in HGHE, an appropriate turbulence model should be selected for calculation. The $k-\epsilon$ equations are the most widely used in numerical computation and engineering applications. They include:

Turbulent kinetic energy equation:

$$\frac{\partial \left(\mathit{\rho k}\right)}{\partial \tau }+\frac{\partial \left(\rho {\mathit{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_i}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\mathit{\rho \epsilon }-Y_M+S_k$$E1

Dissipation rate equation:

$$\frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial \tau }+\frac{\partial \left(\mathit{\rho \epsilon }{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_i}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}+S_{\epsilon }$$E2

where $\rho$ is the fluid density, kg/m³. $\tau$ is the time, s. $k$, $\epsilon$ respectively represents the dissipation rate of turbulent kinetic energy and pulsation energy. $u_i$, $u_j$ is the hourly average speed, m/s. $\mu$, ${\mu }_i$ respectively represent the dynamic viscosity and turbulent viscosity. $G_k$ is the turbulent kinetic energy generation term caused by average velocity gradient. $G_b$ is the turbulent kinetic energy generation term. $Y_M$ is the contribution of pulsation expansion. $C_{1\epsilon }$, $C_{2\epsilon }$, and $C_{3\epsilon }$ are the empirical constants, the values of which are 1.44, 1.92, and 0.09, respectively. ${\sigma }_k$, ${\sigma }_{\epsilon }$ are the Pr number, and the values are 1 and 1.3, respectively. $S_k$, $S_{\epsilon }$ are the user-defined source entries.

The circulating fluid in the buried pipe is water, which is an incompressible fluid and ignores the action of gravity. Thus, it can be simplified as:

Turbulent kinetic energy equation:

$$\frac{\partial \left(\mathit{\rho k}\right)}{\partial t}+\frac{\partial \left(\rho {\mathit{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_i}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k-\mathit{\rho \epsilon }$$E3

Dissipation rate equation:

$$\frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial t}+\frac{\partial \left(\mathit{\rho \epsilon }{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_i}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\frac{\epsilon }{k}\left(C_{1\epsilon }{G}_k-C_{2\epsilon }\mathit{\rho \epsilon }\right)$$E4

Among them, the Prandtl number, turbulent kinetic energy $k$, and dissipation rate $\epsilon$ can be automatically calculated by FLUENT software through fluid physical parameters and velocity.

The standard $k-\epsilon$ model, mass conservation equation, momentum conservation equation, and energy conservation equation together constitute the governing equation of flow and heat transfer, the general form of which is:

$$\frac{\partial \left(\mathit{\rho \varphi }\right)}{\partial \tau }+div\left(\mathit{\rho U\varphi }\right)=div\left({\mathrm{\Gamma }}_{\varphi }grad\varphi \right)+S_{\varphi }$$E5

where $\varphi$ is the universal variable. ${\mathrm{\Gamma }}_{\varphi }$ is the generalized diffusion coefficient. $S_{\varphi }$ is the generalized source term.

The values of parameters in the general equation are shown in Table 2:

Equation	$\varphi$	${\mathrm{\Gamma }}_{\varphi }$	$S_{\varphi }$
Turbulent kinetic energy equation	$k$	$\mu +{\mu }_i/{\sigma }_k$	$G_k-\mathit{\rho \epsilon }$
Dissipation rate equation	$\epsilon$	$\mu +{\mu }_i/{\sigma }_{\epsilon }$	$\epsilon \left(C_{1\epsilon }{G}_k-C_{2\epsilon }\mathit{\rho \epsilon }\right)/k$
Mass conservation equation	1	0	0
X- momentum conservation equation	$u$	$\eta =\mu +{\mu }_i$	$-\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(\eta \frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(\eta \frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial z}\left(\eta \frac{\partial w}{\partial x}\right)$
Y-momentum conservation equation	$v$	$\eta =\mu +{\mu }_i$	$-\frac{\partial p}{\partial y}+\frac{\partial }{\partial x}\left(\eta \frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial y}\left(\eta \frac{\partial v}{\partial y}\right)+\frac{\partial }{\partial z}\left(\eta \frac{\partial w}{\partial y}\right)$
Z- momentum conservation equation	$w$	$\eta =\mu +{\mu }_i$	$-\frac{\partial p}{\partial z}+\frac{\partial }{\partial x}\left(\eta \frac{\partial u}{\partial z}\right)+\frac{\partial }{\partial y}\left(\eta \frac{\partial v}{\partial z}\right)+\frac{\partial }{\partial z}\left(\eta \frac{\partial w}{\partial z}\right)$
Energy conservation equation	$T$	$\mu /\mathit{Pr}+{\mu }_i/\sigma$	0

Table 2.

The value of each parameter in the general equation.

2.3.2 Governing equations of heat transfer of soil

The heat transfer process of soil around the spiral buried pipe also follows the general governing Eq. (5), and the values of its parameters are shown in Table 2.

2.4 Initial and boundary conditions

2.4.1 Upper surface dynamic boundary condition

In this chapter, the data in the Special Meteorological Data Set for the Analysis of the Thermal Environment of Buildings in China are adopted to simulate the hourly meteorological data of the typical meteorological year of Chongqing, as shown in Figure 4, and the dynamic boundary conditions of the upper surface are written into ANSYS_FLUENT software through user-defined function for simulation.

Figure 4.

Hourly meteorological parameters in Chongqing of the typical year.

2.4.2 Initial soil temperature condition

Because HGHE is generally shallow buried depth, shallow soil initial temperature distribution is often uneven due to external influence. Based on this, in order to make the simulation more suitable to the actual project, this chapter selects the soil temperature distribution of the hottest month in Chongqing as the initial temperature condition of the simulation, as shown in Figure 5. Also, the initial soil temperature conditions are written into ANSYS_FLUENT software through user-defined function for simulation calculation.

Figure 5.

Initial soil temperature distribution in the depth direction in the hottest month of Chongqing.

2.4.3 Parameter of soil

The research background of this chapter is in Chongqing, China, and the soil type widely distributed in Chongqing is strongly weathered mudstone, the specific parameters of which are determined by experiment, as shown in Table 3.

Soil type	Strongly weathered mudstone
Density/kg m^−3	2503
Specific heat capacity J kg^−1 °C^−1	1085
Porosity	0.3
Comprehensive thermal conductivity/W m^−1 °C^−1	1.24

Table 3.

Thermophysical parameters of strongly weathered mudstone.

2.5 CFD simulation setup

The HGHE model in this chapter is a three-dimensional model, so the option of 3D model should be selected on the initiator interface. In the calculation of the three-dimensional model, ANSYS_FLUENT provides two different calculation methods, and the two different calculation methods have their own suitable situations. The calculation method of single precision is rough, but the amount of calculation is relatively small, and it is applied to the case that the model is regular. The calculation method of double precision is mainly used for the situation that the length and size of the fluid domain is large, the fluid domain has more than one part, and each domain is connected by the small size pipeline with drastic temperature changes and high thermal conductivity. Therefore, this chapter chooses to use the double-precision solver to solve the problem.

This chapter mainly simulates the flow characteristics of fluid in HGHE, which belongs to the range of incompressible flow fluid, so the uncoupled implicit solution is adopted. Because the heat transfer process in the buried pipe was ongoing all the time and the data fluctuated within the range, the transient computing mode was selected. At the same time, the energy equation option is opened, the $k-\epsilon$ flow model is adopted, and the near-wall function is modified at the same time.

3. Simulation results

The HGHE mainly relies on the temperature difference between the fluid in the pipe and the soil outside the pipe for heat transfer, so the temperature of the fluid in the pipe will directly affect the heat transfer process of the buried pipe, and the temperature of the fluid in the pipe is determined by the inlet water temperature. Therefore, this chapter focuses on studying the influence of different inlet water temperatures on the heat transfer characteristics of the HGHE and the soil thermal response. According to the relevant provisions of the specification Water (Ground) Source Heat Pump Unit (GB/T 19409-2013) and combined with the actual project, four different inlet water temperatures (32, 35, 38, and 41°C) were selected at equal intervals to conduct the simulation research. The spiral tube flow rate was 170 L/h, and the continuous operation was simulated for 120 h. The variation characteristics of soil temperature field and heat flux per unit pipe length with time were analyzed.

3.1 Heat transfer capacity

Figure 6 shows the change curve of heat transfer capacity of HGHE with time under different inlet water temperatures. Figure 6a shows the change of outlet water temperature of HGHE with running time under four different working conditions. It can be found that outlet water temperature gradually increases with time and finally reaches the stable state. Corresponding Figure 6b shows the change of heat flux per unit pipe length with running time under four working conditions. The heat flux per unit pipe length gradually decreases with the passage of time and finally reaches a stable state. Both the outlet temperature curve and the heat flux curve have experienced two stages: the dramatic change stage in the initial stage of operation and the gentle change stage in the later stage. It can be seen that the higher the inlet water temperature is, the more drastic the change at the early stage of operation. This is because the higher the inlet water temperature is, the more heat will accumulate near the spiral pipe in the early stage, which makes the heat exchange resistance of the spiral pipe greater; thus, the heat exchange per unit pipe length decreases more dramatically.

Figure 6.

Heat exchange capacity changing with time of different inlet water temperature conditions. (a) Outlet water temperature changing with time. (b) Heat flux per unit pipe length changing with time.

In order to more clearly explore the influence of inlet water temperature on the heat transfer capacity of HGHE, the relationship curves between heat flux per unit pipe length and inlet water temperature at different times (24, 72, and 120 h) are made, as shown in Figure 7. It can be seen that under the three moments, the heat flux per unit pipe length has a linear increasing trend with the increase of the inlet water temperature. When the heat removal operation reaches 24 h, the heat flux per unit pipe length of HGHE is 5.63, 7.60, 9.56, and 11.52 W/m, respectively, for the inlet water temperature of 32, 35, 38, and 41°C. Under the conditions of four different inlet temperatures (32, 35, 38, and 41°C), the heat flux per unit pipe length of HGHE is 3.83, 5.08, 6.35, and 7.63 W/m, respectively, at the moment of 72 h. When the heat removal operation reaches 120 h, the heat flux per unit pipe length of HGHE is 3.18, 4.22, 5.23, and 6.23 W/m for the inlet water temperature of 32, 35, 38, and 41°C, respectively. It is concluded that when the system runs to the late stable state, the linear increase range of heat flux per unit pipe length of HGHE becomes small with the inlet water temperature. For example, at the moment of 120 h, the heat flux per unit pipe length increases only about 0.34W/m when the inlet temperature increases by 1°C.

Figure 7.

Heat flux per unit pipe length changing with inlet water temperature.

3.2 Thermal response of soil

Figure 8 shows the temperature distribution field of half longitudinal sections under different inlet water temperature conditions. It can be seen from the figure that the higher the inlet water temperature is, the more obvious the soil heat accumulation phenomenon near the spiral buried pipe is, and the stronger the thermal interference effect is. In particular, the temperature of the soil inside the HGHE increases dramatically with the increase of the inlet water temperature, which is mainly due to the limited volume of soil in the pile. With the progress of the heat removal condition, the heat discharged from the spiral pipe around HGHE will be transferred to the pile, resulting in heat accumulation, while the heat inside HGHE can only be discharged through the upper and lower opening surface. Therefore, the increase of inlet water temperature means that more heat will accumulate inside HGHE, so the thermal interference inside the pile is more and more serious, which will also greatly hinder the improvement of the heat transfer capacity of HGHE. It can also be seen that the increase of inlet water temperature has little influence on the thermal action range of HGHE, but mainly increases the thermal response temperature of soil within the thermal action range. Therefore, it shows that increasing the inlet water temperature of HGHE has no significant effect on improving the heat exchange performance. On the contrary, it will increase the outlet temperature of HGHE, resulting in higher condensing temperature and lower the efficiency of the GSHP system.

Figure 8.

Soil temperature distribution of different inlet water temperature conditions. (a) 32°C, (b) 35°C, (c) 38°C, (d) 41°C.

4. Model validation

In order to further verify the accuracy of the numerical model established in this chapter, a heat transfer experimental platform corresponding to the HGHE was built, as shown in Figure 9. The verification experiment was carried out, and the measured project was compared with the simulated results, through which the verification of the model could be completed.

Figure 9.

Heat transfer experimental system of HGHE.

As can be seen from the comparison curve of inlet and outlet water temperature in Figure 10, the numerical simulation results are close to the experimental results, and the maximum difference is less than the 5% relative error commonly used in engineering. But the simulation results are slightly larger than the test results. In practice, there may be a series of conditions such as underground water migration, which are simplified in the model, so there are differences, but the differences are not large, resulting in the measured outlet temperature that is slightly lower than the simulated outlet temperature. In general, the operating results of the model are close to the measured values of the actual project, and the model is basically consistent with the actual values, indicating that the model is basically correct and has a certain guiding role for the actual project.

Figure 10.

Comparison of simulation and experiment results.

5. Conclusions

This chapter focuses on the application of CFD technology in GSHP system. Firstly, for the simulation research field of GHE heat transfer characteristics and heat transfer process, the research progress of scholars from various countries is reviewed. Then, in order to further elaborate and introduce the application of CFD technology in GSHP system, the HGHE is taken as an example. Considering the dynamic environment conditions and the initial soil temperature conditions, the ANSYS_FLUENT software commonly used in CFD technology was used to simulate and study the heat transfer characteristics of HGHE under different inlet water temperature conditions, and the accuracy of the numerical model was verified by experimental means. The following conclusions are reached through the research:

1. CFD technology can accurately reflect the actual heat transfer process of HGHE, and it has a certain guiding role for practical engineering.

2. The heat flux per unit pipe length increases linearly with the increase of the inlet water temperature. When the inlet water temperature of HGHE is low, the heat exchange of the buried pipe can be increased by increasing the inlet water temperature. However, when the inlet water temperature of HGHE is high, continuing to increase the inlet water temperature will seriously aggravate the soil thermal interference near the spiral pipe, which will greatly hinder the heat transfer process of HGHE. At the same time, too high inlet water temperature will increase the outlet water temperature, indirectly causing higher condensation temperature, resulting in lower efficiency of the GSHP system. Therefore, the thermal interference effect and the condensing temperature of GSHP system should be comprehensively considered in the actual project, and the appropriate inlet water temperature of HGHE should be selected.

Acknowledgments

This work received support from Natural Science Foundation of China [52106283].

Conflict of interest

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work. There is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, our work.

Appendices and nomenclature

T	temperature (°C)
x, y, z	three-dimensional direction (m)
ui, uj	hourly average speed (m s−1)
μ, μi	dynamic viscosity and turbulent viscosity
Gk	turbulent kinetic energy generation term caused by average velocity gradient
Gb	turbulent kinetic energy generation term
YM	contribution of pulsation expansion
C1ε, C2ε, C3ε	empirical constant, the values of which are 1.44, 1.92 and 0.09 respectively
Sk, Sε	user-defined source entries
ϕ	universal variable
Γϕ	generalized diffusion coefficient
Sϕ	generalized source term
Greek symbols
ρ	density (kg m−3)
τ	time (s)
k, ε	dissipation rate of turbulent kinetic energy and pulsation energy
σk, σε	Pr number, and the values are 1 and 1.3 respectively