Heat and Mass Transfer Behavior Prediction and Thermal Performance Analysis of Earth-to-Air Heat Exchanger by Finite Volume Method

: A comprehensive numerical study on coupled heat and mass transfer in an earth-to-air heat exchanger (EAHE) is conducted by self-complied program based on the ﬁnite volume method. The soil thermal and moisture coupled characteristics in the vicinity of the pipe and the thermal performance of the EAHE are evaluated by a two-dimensional simulation model. The model of the EAHE is veriﬁed by the experimental data, which achieved a good agreement with each other. The numerical results show that there is an obvious moisture peak in the radial direction, and the peak position radially moves away from the wall of the pipe over time. It is also found that the thermal performance of the heat and mass transfer model in soil is better than the pure heat conduction model.


Introduction
Nowadays, the energy consumption is growing rapidly. Furthermore, in last decades, a large amount of energy consumption comes from non-renewable fossil fuels such as coal, oil, and natural gas. With the increasingly global energy crisis and the deterioration of the environment, developing and exploiting new green energy have been received worldwide attention. Obviously, solar energy is exactly the ideal type which is inexhaustible and pollution-free. A large amount of solar radiation is stored at superficial layers of the ground and it decreases rapidly with the increasing depth of the soil. As a result, the soil temperature under 4 m is nearly constant [1]. Utilizing the geothermal as a heat source and sink, the EAHE can heat and cool the ambient air for reducing the energy consumption of buildings. The use of EAHE for cooling and heating agricultural greenhouses and buildings has gained much achievement in the past several years [2,3].
A widely used physical model for studying coupled heat and mass transfer in a porous medium was presented by Philip [4] and De Vries [5,6]. In a porous material, the differential equations for the coupled heat and mass transfer process have been based on temperature and moisture concentration gradients as driving potentials, which means that the moisture movement in the porous material is driven by the temperature gradient. Simultaneously, the moisture will trend to redistribute itself under the created moisture gradient. The total moisture flux was divided into two components: liquid water and water vapor. An alternative model was given by Eckert et al. [7,8] in 1979. Jahangir [9] developed the new coupled differential equations for heat, mass and air transfer process in porous medium based on the above alternative model, which could result in a better accuracy and significantly reduce the simulation time by concentrating the computational effort in the vicinity of the soil-pipe.
Shah et al. [10] found that liquid moisture moved from the cold end to the warm end while water vapor migrated in the opposite direction. They concluded that the primary cause of liquid water migration was the moisture concentration gradient, while the water vapor migration was predominantly due to the difference in the soil temperature.
Through the experiment in summer, Balghouthi et al. [11] studied the thermal and moisture characteristics of dry and wet soil which are heated by buried capillary plaits under three different operational conditions: heating at 70 • C, heating at 40 • C and without heating. It was concluded that the surface temperature amplitude in moist soil was superior to dry soil. In order to determine the thermal and mass diffusion coefficients for two different types of soils, Shah et al. [10] conducted a 1-D steady state soil column experiment. They found that the ratio of the thermal diffusion coefficient to the moisture diffusion coefficient monotonically increases with increasing the moisture concentration, and then gradually decrease after reaching a peak value. To measure the primary and cross-coupling transport coefficients, Jury [12] conducted a two-column experiment through sand for the coupled and simultaneous heat and moisture transfer process.
In order to determine the optimal design of a closed loop EAHE for greenhouse heating, an exergo-economics analysis was conducted by Ozgener et al. [13][14][15][16][17]. Thanu et al. [18] conducted an experimental study of an EAHE system in single pass mode at Guomohar farmhouse. They established the correlations between outlet air and inlet air temperature, indicating the high effectiveness of the EAHE system in summer. With the application of an artificial neural network concept Kumar et al. [19] proposed two models, deterministic and intelligent, that can help the designers evaluate the heating and cooling potential of an EAHE. Furthermore, they found that the deterministic model with an accuracy of 75.3%, while the intelligent model accuracy was 72.6%. Kumar et al. [20] presented a transient model that can evaluate thermal performance and dynamic performance of building and the energy conservation potential of an EAHE by using MATLAB (R2015b, MathWorks, Natick, MA, USA) and the finite difference method. In [21], the thermal performance of the EAHE was experimentally investigated in Brazil. They found that such a system could achieve about 48% of cooling or heating demands. A transient analytical model of an EAHE connecting to a wind tower was presented by Benhammou et al. [22]. The results showed that the thermal performance of the wind tower coupled to EAHE was superior than the conventional cooling tower. In the article [23], Kaushal et al. utilized a commercial CFD to predict the thermal performance of hybrid EAHE and RSM (response surface methodology) to optimize the process parameters. The results revealed that the thermal performance of hybrid EAHE was better than individual EAHE. To evaluate the thermal performance of EAHE under transient operational conditions, Bansal et al. [24] proposed a new concept of derating factor, which is determined by the EAHE air temperature drops in both transient condition and steady state condition. The results indicated that, in transient state, the thermal performance of EAHE gradually deteriorates under the continuous operation.
To present an investigation of the soil thermal and moisture coupled characteristics in the vicinity of the pipe under warming conditions and the thermal performance of the EAHE in transient condition, a transient bi-dimensional model is established in this paper. In addition, a comparative research is carried out for investigating the influence of different heat transfer models, which focuses on the thermal performance parameters of EAHE as follows: the daily mean efficiency, the difference between suction air temperature and delivery air temperature, the daily cooling potential, and derating factor in transient condition.

Description of the System
The major component of the EAHE (Figure 1) is a buried PVC (poly vinyl chloride) pipe in the soil. The length and the interior diameter of the buried pipe is L and D, respectively. As shown in the figure, during the operation of the system, the outlet of the pipe is connected to an agricultural greenhouse while the inlet is combined with ambient air. At an adequate depth, the soil temperature is comparatively constant during the whole year due to the delayed thermal response of soil, and it is higher than the ambient air temperature in winter but lower in summer. Accordingly, passing through the buried pipes, the ambient air will be heated in winter and cooled in summer. After being injected into the agricultural greenhouse and mixed with the indoor air, the ambient air can make the air temperature and humidity of the agricultural greenhouse meet the basic needs of the crop growth.

Description of the System
A widely used formulation for heat and mass transfer in porous material is written as follow, in which the key assumptions are adopted [25]: 1.
The total system pressure is approximated by air pressure 2.
The vapor-liquid interface is a function of temperature only 3.
The total pressure is constant during the transfer process 4.
Heat transfer by radiation is negligible 5.
The air is incompressible and the property of air is constant The conservation equations for energy and mass within the soil in 2D cylindrical coordinates can be stated as follows: where the coefficient of the above equations are summarized in Table 1.
Thermal liquid diffusivity D T,w m 2 /(s K) Thermal vapor Diffusivity D T,vap m 2 /(s K) The governing transport equations for air in the pipe in 2D cylindrical coordinates are given below: Momentum equations: Energy equation: Transport equations for the k-ε model: The boundary conditions of the computational domain are: Furthermore, the outlet air velocity in the axial direction must satisfy total mass conservation laws. The initial conditions of the computational domain are:

General Remarks
The aforementioned conservation equations are discretized by the finite volume method on a staggered grid system [26,27]. The QUICK (Quadratic Upstream Interpolation for Convection Kinetics) scheme and central difference scheme are used to discretize the convection term and the diffusion term, respectively, and an efficient segregated algorithm, SIMPLER (Semi-Implicit Method for Pressure Linked Equations Revised), is applied to the pressure-velocity solution. The algebraic equations for the whole computational domain can be solved by the TDMA (tridiagonal matrix algorithm) +ADI (alternating direction implicit) method. The linearization method is used to deal with the source terms of the governing equations above. The full-field computation method is adopted to solve the fluid-solid coupling heat transfer process. Accordingly, the dependent variables of soil and air region are solved simultaneously. In order to guarantee the continuity of the heat flux at the soil-pipe interface, ρc p is regarded as a general density in the energy equation [28]. In the soil region, by adopting a very large value of physical viscosity in the momentum equations, for example 10 30 , the zero velocity of soil area can be gained. To ensure the convergence of iteration, under-relaxation of the dependent variables and pressure can be incorporated into the solution process of the algebraic equations. According to the upwards treating processes, a 2-D program has been developed to solve the model.

The Discretization of Governing Equation
A general expression for the afore mentioned governing equations can be written as: where Γ Φ is the nominal diffusion coefficient and S Φ is the source for the general variable Φ. Equation (11) is discretized by the finite volume method on a staggered grid system as shown in Figure 2. To enhance the convergence rate, the source term S Φ is linearized as follows: where S P ≤ 0, which guarantees the requirement of diagonal dominant. For two-dimensional problems, after integration over the control volume, the discretized equation can be obtained as follows: The variables Φ n,e,w,s at the control interfaces can be determined using an interpolation or extrapolation involving the values of the neighboring grid points. F n,e,w,s are the flow rate at the interfaces. D n,e,w,s are diffusive conductance. ∆V is the volume of the control volume. Using Patankar's notation [26], Equation (13) can be written in the form: where a P , a 0 P and a N,W,E,S are coefficients of the discretized equation. b is constant term in the discretized equation.
Underrelaxation of the dependent variable is incorporated into the solution process of the algebraic equations, Equation (14) becomes: where χ is the underrelaxation factor. Φ 0 P is the value of general variable Φ of the previous time level.

Treatments of Near Wall Nodes
The correct processing method of the near wall nodes is an important step for the programming calculation. In this paper, the wall function method is applied to establish the fluid-solid interface boundary conditions. As shown in Figure 3, some variables of the first inner node P can be stated as follows: 1.
The equivalent viscosity conductivity between the first inner node and the wall is, The equivalent thermal conductivity between the first node and the wall is,

3.
By adopting large coefficient method, the dissipation rate of the first inner P is taken the given Ky P

4.
The first order normal derivatives of velocity normal to wall and the turbulence kinetic energy equal to zero.
where y P + = ρy P c µ

Cooling Potential
The cooling potential of EAHE for a period of time τ can be calculated by the following equation: T a (in) (j) and T a (out) (j) are the air temperature of each grid point at inlet and outlet surface of the pipe respectively. u in (j) and u out (j) denote the air velocity of each grid point at inlet and outlet surface of the pipe, respectively.

Average Nusselt Number of the Wall
The equivalent average Nusselt number of the wall can be calculated according to the following correlations: where The temperature of n-interface in Figure 2 can be calculated by the interpolation from neighboring grid points according to the heat flux equilibrium.
T sw is the average soil temperature at the soil-pipe interface.
S sur f ace is the area of the soil-pipe interface. The R N and R P are radial distance of N point and P point, respectively.

Mean Efficiency
The mean efficiency of EAHE for a Period τ is calculated by the following relation, are the average air temperature at inlet and outlet surface of the pipe respectively. T s is soil temperature.

Derating Factor
Temperature drops gained by EAHE in steady state condition and transient condition are used to determine the derating factor which is defined as follows: DF l,t : derating factor at a distance of 'l' from EAHE inlet, after time 't'. T a (l) is the air temperature at length of 'l' m from the pipe inlet.

Converge Criteria
In our previously numerical simulation about heat and moisture transfer in soil by looking into the process of convergence, it has been verified that when the two convergence criterion as follows are satisfied at the same time, the convergent solution can be well obtained: m ≤ 10 −6 , and, N u,ave q+200 − N u,ave q /N u,ave q+200 ≤ 10 −6 (25)

The Structure of Computation Program
The numerical simulation program of heat and moisture transfer in soil is divided into two parts: main program and user subroutine. The main program determines the flow of the entire numerical simulation and the operation and termination of the program. The calculation conditions and contents are completed by different subroutine modules. The flowchart of the computation program is shown in

Grid-Independence Test
In this paper, five sizes of grid were established to determine the accuracy of the average Nusselt number of the wall (32 × 90, 64 × 180, 128 × 360, 192 × 540, and 256 × 720), and the results are shown in Figure 5. The value of the average Nusselt number of the wall changes no more than 2%, when the total grid number increases from 4.6 × 10 4 to 1.8 × 10 5 , which indicates a great stability of the grid. In order to save computation time and guarantee the accuracy of the simulation, the optimum grid size, 128 × 360, is adopted in the following simulation.

Self-Complied Program Examination
The accuracy and feasibility of self-complied program was verified through experimental data. A long-term experimental test was carried out by Liu [29] from 25 April 2014 to 27 August 2014 using the existing facility which contains an 18 m heat exchanger pipe, the diameter of which is 110 mm, and the pipe is buried under the ground of the agricultural greenhouse for the depth of 1 m. The inlet velocity of the air was 4.5 m/s. The temperature of soil at different depth were measured. In addition, the schematic diagram of this experimental platform is shown in Figure 6. In order to clearly compare the measured data with simulated results concluded from the self-complied program, a whole day's experimental data of 8 June 2014 was adopted to verify the model. The comparative results are presented in Figures 7a and 8a, and the error analysis of model estimation was shown in Figures 7b and 8b. Figure 7a illustrates the variation of the measured outdoor air temperature, outlet air temperature and simulated outlet air temperature with time, indicating the simulated outlet air temperature is in good agreement with the experimental data. Furthermore, Figure 7b demonstrates the matching level of testing and simulated data, which indicates the mean square error (MSE) is 0.73 • C and the mean relative error (MRE) is 8.5%. The comparison of the soil temperature along the radial direction between the testing and simulated data is shown in Figure 8a, which indicates that the two sets of data match very well. As shown in Figure 8b, the MSE is only 0.57 • C and the MRE is 6.4%. The error can satisfy the needs of the engineering practice. According to the comparison and error analysis of testing and simulated data, the self-complied program used in this paper can reliably and precisely deal with the 2D coupled heat and mass transfer in soil.

Soil Temperature and Moisture Concentration Analysis in Radial and Axial Directions
The soil temperature and moisture concentration profiles from 0 m to 0.5 m away from the soil-pipe interface are simulated and plotted in Figure 9a. The initial soil moisture concentration is assumed to be 20%, and the air flow rate and inlet air temperature are set as 154 m 3 /h and 38 • C, respectively. As shown in Figure 9a, after 24 h continuous operation, the soil moisture concentration monotonically first increases, and then gradually drops to initial state after reaching the moisture peak value. As a result, the moisture concentration curve at the end of 24 h can be separated into three regions: temperature driven section, common driven section, and undisturbed section.
In the first region, moisture migrates under temperature gradient from the higher temperature regions to lower temperature regions, while simultaneously tending to redistribute itself in reverse under the created soil moisture concentration gradient. The main cause of soil moisture migration is temperature gradient in this region. Accordingly, the soil moisture concentration drops from 20% to 18.81% at the soil-pipe interface and has been below initial value in the first 0.13 m of the first region. The soil moisture concentration increases monotonically from the soil-pipe interface and exceeds the initial value about 20% at the last part of the first region due to the accumulation of moisture concentration. Finally, the soil moisture concentration peaks at 20.16% when the temperature gradient and the created moisture gradient reach a temporary equilibrium at 0.18m away from the soil-pipe interface.
In the second region, soil moisture concentration gradually decreases from the peak value to the initial value. Differing from the first region, the moisture migration direction caused by the temperature gradient is same as the moisture migration direction caused by the moisture gradient. The moisture radially migrates outwards the soil-pipe interface under both moisture gradient and temperature gradient. Accordingly, this region is called common drive section.
The third region is radially away from the soil-pipe interface from 0.37 m to 0.5 m. In this undisturbed region, moisture concentration remains nearly constant. From Figure 7a it can be seen that soil moisture concentration reaches the peak value of 20.07% at the radial distance of 0.13 m after 12 h of continuous operation. As a result, it is concluded that the peak value increases and the peak position moves radially outwards with time, and the area of the undisturbed region decreases gradually with time.
Although the soil temperature decreases gradually along the radial direction from the soil-pipe interface, the decreasing rate drops with the increase of the distance away from the soil-pipe interface. Noticing the soil temperature increases with the continuous operation of the system. From Figure 9a, the penetration of heat at the entrance region is not beyond four times of the pipe diameter at the end of 24 h. As a result, the computational domain extends to five times of the pipe diameter in this paper. Figure 9b represents the change of the soil temperature and moisture concentration away from the soil-pipe interface about 0.05 m along with the pipe length. It is easy to find that the soil temperature gradually decreases along the length of the pipe, and the moisture concentration increases rapidly near the air entrance region and then gradually slows down. Meanwhile, the soil moisture concentration becomes smaller with time.

Thermal Performance Evaluation under Different Heat Transfer Models
The effect of the pipe length on the daily mean efficiency is illustrated in the Figure 10. It is found that the daily mean efficiency increases with the increasing of the pipe length for different heat transfer models. For a change in the pipe length from 5 to 15 m, there is an obvious increase in daily mean efficiency of heat and mass transfer model from 34.2% to 72.5% and a variation gained from pure heat conduction model from 28.3% to 62.8%. As a result, it has been proved that the heat and mass transfer model can yield the best daily mean efficiency is compared to the pure heat conduction model.  Figure 11 shows the evolution of difference between temperature of inlet air and outlet air as function of operating duration. The value of air temperature difference between inlet and outlet varies between 15.3 • C and 13.8 • C for heat and mass transfer model and pure heat conduction model after 12 h of continuous operation. The value of the same is 14.4 • C and 12.2 • C when EAHE runs continuously for 24 h. According to the graph, the air temperature difference between inlet and outlet drops for the long period of continuous operation. However, the temperature difference between two heat transfer models becomes bigger with time. This can be explained by the fact that moisture migration leads to the soil temperature in the vicinity of buried pipe of heat and mass transfer model is lower than the pure heat conduction model during the running process of the system. Accordingly, the outlet air temperature of heat and mass transfer model is lower than the pure heat conduction model. Figure 12 represents the variation of the daily cooling potential along with the pipe length for different heat transfer models. It is shown that the daily cooling potential increases as function of pipe length for two heat transfer models. For which, it is observed a rise from 8.4 to 18.1 kW h corresponds to variation of pipe length from 5 to 15 m for heat and mass transfer model. For the same change in the pipe length, it is observed a rise from 7.2 to 15.9 kW h for pure heat conduction model. And, the difference of the daily cooling potential of two heat transfer models becomes bigger according to the pipe length. Additionally, it is true that increasing the pipe length has a positive effect of increasing the daily cooling potential, but the increasing rate drops according to the pipe. Figure 11. The difference between inlet air temperature and outlet air temperature for different heat transfer models. In Figure 13, the evolution of the derating factors at different lengths of pipe after different times of continuous operation under transient conditions for two heat transfer models is shown. It is mentioned that the derating factor decreases as function of pipe length. The derating factor drops from 0.22 to 0.1 and 0.72-0.31 after 1 and 24 h of continuous operation of EAHE, respectively, for pure heat conduction model. The value of the same for heat and mass transfer model is 0.06-0.16 and 0.25-0.58. The value of the derating factor at a distance of 6 m from the pipe inlet for continuous operating duration of 24 h varies between 0.11 and 0.43, 0.13 and 0.52 for the heat and mass transfer model and the pure heat conduction model, respectively. Therefore, it is concluded that thermal performance deteriorates with the duration of continuous operation. Larger values of the derating factor mean greater deterioration in the thermal performance of EAHE under transient conditions. As shown in Figure 13, the value of the derating factor of pure heat conduction model is greater than heat and mass transfer model for an equal period of time. Furthermore, curves of the derating factor for pure heat conduction model are much steeper than the corresponding curves for heat and mass transfer model. According to the above analysis, deterioration in thermal performance of EAHE for pure heat conduction model is larger and more rapid than the heat and mass transfer model.

Conclusions
In this paper, a numerical study conducted on an EAHE integrated with an agricultural greenhouse is presented. The thermal performance of EAHE under different heat transfer models and soil temperature and moisture behaviors are discussed in detail. Some main conclusions are summarized as follow: 1.
There is an obvious moisture peak in the radial direction, which has an important effect on the heat and mass transfer process. Furthermore, the location of peak will move radially outwards and the peak value becomes bigger with time.

2.
Compared with the pure heat conduction model, the lower outlet air temperature and soil temperature can be obtained for the heat and mass transfer model for an equal period of continuous operation duration.

3.
Heat and mass transfer model can yield best mean efficiency, daily cooling potential and thermal performance of EAHE in transient conditions are compered to pure heat conduction model.

4.
Thermal performance of EAHE deteriorates with duration of operation. Accordingly, understanding the soil thermal saturation and recovery of an EAHE system under different types of operation (intermittent or continuous) is very important.