An improvement to the simulation of phase-change heat-transfer during soil freezing and thawing

M INING SCIENCE AND TECHNOLOGY Mining Science and Technology 19 (2009) 0262–0268

www.elsevier.com/locate/jcumt

An improvement to the simulation of phase-change heat-transfer during soil freezing and thawing YANG Tao, ZHENG Mao-yu School of Municipal and Environment, Harbin Institute of Technology, Harbin, Heilongjiang 150090, China Abstract: A new variable time step method, which is called the backwards calculating time step method, is presented in this paper. It allows numerical simulation of soil freezing and thawing while avoiding “phase change missing and overflowing”. A sensitive heat capacity model is introduced through which the calculation errors are analyzed. Then the equation using the self-adjusted time step is presented and solved using finite differences. Through this equation, the time needed for a space cell to reach the phase change point temperature is calculated. Using this time allows the time step to be adjusted so that errors caused by “phase change missing and overflowing” are successfully eliminated. Above all, the obvious features of this method are an accelerated rate for adjusting the time step and simplifing the computations. An actual example proves that this method can accurately calculate the temperature fields during soil freezing and thawing. It is an improvement over traditional methods and can be widely used on complicated multi-dimensional phase change problems. Keywords: backwards calculating time step method; phase change missing; phase change overflowing

1

Introduction

The development of science and technology has caused phase change heat transfer problems to arise in the many domains. For example, the freezing method is a prime construction technology for mine-building in deep alluvium. It is very important to develop numerical simulations of the temperature fields during soil freezing and thawing to allow accurate estimations of the frozen wall[1–2]. The freezing and thawing of soil is a heat conduction problem involving both solid and liquid phases and is an example of the well-known Stefan problem. Since the process is very non-linear it is impossible to solve analytically so it can only be estimated by numerical calculation. The common numerical calculation methods used for heat transfer are mainly: the finite difference method, the finite element method and the finite boundary element method[3]. Although the discretization used in these various methods is different all replace the continuous physical measurement in time and space coordinates with limited, discretized values at space nodes, along with algebraic equations defined at the nodes. If the space coordinates are constant and the time step is selected inappropriately, “phase change missing and overflowing” can easily occur.

The commonly used numerical calculations for solving the temperature fields during soil freezing and thawing are the variable space-step method and the variable time-step method[4–6]. The former considers a grid transformation or trace of the moving interface between grid node coordinates and is more complex. The latter reduces the size of the time step as the soil temperature approaches the temperature where phase change occurs. Although simple, the latter method does not accurately adjust cell temperatures to the phase change point and, furthermore, computing times for the method are too long. Reference [7] describes a method in which completely and partly implicit difference schemes were alternated. This method needed to adopt different discretization approaches. The soil-phase-increasing method assumed a phase change increment and calculations were repeated until the phase change increment reached the true value[8]. This method was unable to compute the time of phase change, and computation times were too long. Reference [9] shows that different temperature ranges for the phase change would have an effect on the final results. Phase change errors were not sufficiently reduced by setting a lager temperature range for the phase change. As a result, a simple, convergent and precise calculation method is needed to solve these complex nonlinear problems.

Received 24 June 2008; accepted 11 September 2008 Project 2006G1662-00 supported by the Key Science and Technology Project of Heilongjiang Province Corresponding author. Tel: +86-13796632173; E-mail address: [email protected]

YANG Tao et al

An improvement to the simulation of phase-change heat-transfer …

A self-adjusted time step method, called the backwards calculating time step method (BCTM) is presented here. It is based on finite differences and the sensitive heat capacity mode. This method, in contrast to the traditional variable time step method, can successfully avoid “phase change missing and overflowing” and enhance the rate at which the time step is adjusted. The method also achieves higher precision. A two-dimensional model of heat transfer was established to simulate the freezing and thawing of the soil. Since the simulation results of the BCTM are consistent with the test results the method has been proven feasible.

2 Mathematical model The common multi-dimensional unsteady-state differential equation describing heat conduction during the freezing and thawing of soil is

C

∂θ ∂T = ∇ [ λ∇ T ] + L ρ s s ∂t ∂t

(1)

temperature, t is time and L is the freezing latent heat of the soil water. During freezing, the unfrozen water content is a function of the temperature

θ u = θ u (T )

(C + L ρ w

C*

­λf ° λ − λu ° * λ = ®λf + f ªT − (Tf − ΔT ) º¼ 2ΔT ¬ ° °¯λu

Phase change missing and overflowing

The precision of common finite difference schemes falls in the order: CN more precise than implicit, which is more precise than the sensitive scheme. However, the implicit scheme’s stability is the best[10].

(3)

∂T = ∇ ª¬λ *∇T º¼ ∂t

­Cf ° Lρ φ C + C ° u C* = ® w + f Δ 2 T 2 ° °¯Cu

C and λ are the volume thermal capacity of the soil and the thermal conductivity of soil, T is the transient

3

∂θ u ∂T ) = ∇ [ λ ∇T ] ∂T ∂t

The derivative term relates the latent heat to time and brings difficulties for the solving of the temperature fields. The phase change of soil water does not strictly occur at a specific temperature, but occurs over a small temperature range. Therefore, the latent heat of soil water can be considered to be a sensitive heat capacity over this small temperature range. The above problem is transformed into the “single-phase” nonlinear heat conduction problem and thus Eq.(3) can be written as

the soil water content by volume, soil unfrozen water content and soil ice content, respectively. In addition, ρ w and ρs are soil water density and ice density,

Cu are volume thermal capacities of frozen and unfrozen soil, λf and λu are thermal conductivities of frozen and unfrozen soil, φ is the soil water content by volume and Tf − ΔT ≤ T ≤ Tf + ΔT is an assumed small temperature range for the phase change. Temperature is divided into three parts: a non-phase change range, a rapid phase change range, and a completed phase change range. Tf − ΔT and Tf + ΔT are the phase change points.

(2)

θs in Eq.(1) is eliminated, and the derivative rule for a compound function is used, to give

ρ where θ = θ u + s θ s , and θ , θ u and θs are ρw

where C * is the equivalent volume thermal capacity, λ * is the equivalent thermal conductivity, Cf and

263

(4)

where (T < Tf − ΔT ) (Tf − ΔT ≤ T ≤ Tf + ΔT ) , (T > Tf + ΔT )

(T < Tf − ΔT ) (Tf − ΔT ≤ T ≤ Tf + ΔT ) (T > Tf + ΔT )

Eq.(4) is solved by the implicit scheme. Through in* putting the initial temperatures and assuming λ and C * at the next time increment, the temperatures at the next time increment are given. Soil is a multi-phase material, which mainly contains water and solid minerals. When the soil temperature drops to the phase change point the temperature at the cell changes from T 1 > Tf + ΔT to T 2 < Tf + ΔT within the time step. The values of λ and C

*

*

vary within the same step. If T 1 is used to *

determine λ * and C * , the cell has a C that neglects the phase change latent heat term: This is “phase change missing.” Also, λ * is wrong. If T 2 is used to determine λ * and C * then the cell has a C * that over-estimates the phase change latent heat term: This is “phase change overflowing” again, λ * is also wrong. The error in C * is called “phase change missing and overflowing.” The value of this

264

Mining Science and Technology

α PT k − (α ET k + α WT k + α NT k + α ST k ) P

=α T

where α E =

αN =

(5)

(6)

rp Δr

( δ x )e λ

k −1 e

, αW =

rp Δr

(δ x )w λwk −1

,

rs Δx rn Δx , αS = , k −1 (δ r )n λn (δ r )s λsk −1

(α p − α 0 )T − (α T + α T + α T + α T ) − S 2 − S3 k E E

S

+ S 2 + S3

Tpk −1 − Tpk k p

N

α0 =

Eq.(5) is discretized by the method of completely implicit finite differences. Energy of absorption or

Δt =

W

C *k −1ΔV , α P = α E + α W + α N + αS + α 0 Δt where the superscripts k and k–1 denote kƸt and (k–1)Ƹt moments in time, respectively, and the subscripts E, W, N and S denote cells adjoining cell P. S2 and S3 are the rates of energy charge or discharge under the second and third boundary conditions, respectively. For discretizations performed under different coordinate systems refer to Reference [11]. The time equation for cell P as its temperature k −1 to Tp is derived from Eq.(6) as passes from Tp

The steps involved in BCTM are: estimating if “phase change missing and overflowing” will occur; calculating the time needed for a cell temperature to reach the phase change point; adjusting the time step accordingly; and, finally, forcing the cell temperature to fall on the phase change point after re-iteration. Specifics follow. Eq.(5) describes dynamic heat conduction in cylindrical coordinates

· 1 ∂ § * ∂T · ¸ + ⋅ ¨ rλ ¸ ∂r ¹ ¹ r ∂r ©

E

k −1 0 P

4 Backward calculating the time step

∂T ∂ § * ∂T = ¨λ ∂t ∂x © ∂x

No.2

loss under the second and third boundary conditions is considered to be the equivalent of a heat source in the control cell close to the boundary[11]. Then Eq.(5) can be written as

error can be called the phase change error. During numerical calculation the cell temperature belongs to different temperature ranges before and after an iteration step. How λ * and C * change can not be estimated only through cell temperature. The phase change latent heat term must be included: without it the appearance of a phase change error is inevitable. If the error is ignored, the calculated temperature will either be too small or too large. To overcome these difficulties, the BCTM is presented in this paper.

C*

Vol.19

k W W

Suppose that phase change occurs in cell P when the time changes from k–1 to k, and that Tm

k N N

k S S

C * ΔV

(7)

The time equation for cell P as its temperature passes from Tpk −1 to Tm should then be

( Tpk < Tm < Tpk −1 ) is the phase change temperature.

Δt ′ =

Tpk −1 − Tm (α p − α 0 )Tm − (α T + α T + α T + α T ) − S 2 − S3 k E E

k W W

where Δt ′ is the time step calculated backward. Eq.(6) is now re-calculated at the k–1 moment. Δt ′ is iterative time and δ is calculation error. While the soil is freezing and Tm − δ < Tp ≤ Tm , or while

k N N

k S S

C * ΔV

(8)

tion turns to calculating the next phase change interface. The time equation of cell P as the temperature k −2 passes from Tp to Tpk −1 as found through Eq.(6) is

the soil is thawing and Tm ≤ Tp < Tm + δ , the simula-

Δt =

Tpk − 2 − Tpk −1 (α p − α 0 )T

k −1 p

− (α T

k −1 E E

+α T

k −1 W W

+α T

k −1 N N

+α T

k −1 S S

) − S 2 − S3

C *ΔV

(9)

Assuming

Tpk − 2 − Tpk −1 ≥ Tpk −1 − Tm

(10)

The time equation for cell P as its temperature passes from Tpk −1 to Tm is

Δt ′ =

Tpk −1 − Tm (α p − α 0 )Tpk −1 − (α ETEk −1 + α WTWk −1 + α NTNk −1 + α STSk −1 ) − S 2 − S3

C *ΔV

(11)

YANG Tao et al

An improvement to the simulation of phase-change heat-transfer …

Eq.(11) must satisfy the Eq.(10). The closer k −1 (Tpk − 2 − Tm ) / 2 is to Tp , the faster the convergence of Eq.(11) will be. When calculating soil freezing and thawing, or while investigating other material phase changes, the largest time step is selected from to the two time steps calculated by Eqs.(11) and (8). The numerical analysis steps are as follows: 1) Judging Suppose that during iteration the soil temperature gradually drops until soil begins freezing. Subscripts i and j denote two minimum temperature cells at mok −1 ments k–1 and k. At those times, Ti ,min and T jk,min denote the temperature of cells i and j. The initial it* eration temperature, T1, is used to determine λ and k −1 C * . When Ti ,min > Tf + ΔT and T jk,min < Tf + ΔT “phase

change

k −1 Tf − ΔT ≤ Ti ,min

missing” and ≤ Tf + ΔT

exists.

When

T jk,min < Tf − ΔT

“phase change overflowing” exists. It is certainly possible that parts of cells appear “Phase change missing” while other parts of cells simultaneously appear “Phase change overflowing”. 2) Computing Δti Suppose cell i appears to be “Phase change missing” at time step k. Δti is the time needed for cell i, whose temperature changes, to reach the phase change point from the k-1 step. The time step is modified by Eq.(8) or (11), and a method of computk −1 ing the time stage by stage is adopted, to allow Ti ,min to approach Tm ( Tm = Tf + ΔT ) gradually. To make the temperature of cell i lie in the phase change range at the beginning of the next iteration, the error δ is set so Tm = Tf + ΔT − δ . The method is: In the first stage, the times Δt1′ and Δt1′′ of cell i as the temperature goes from Tpk −1 to Tm are calculated

from Eqs.(8) and (11), respectively. The larger one is selected as Δt1 , and then Δt1 is put into Eq.(6) as the iterative time at the k–1 step. The cell i temperak −1 to Ti1 . ture is then computed to change from Ti ,min If Ti1 < Tf + ΔT − δ , the smaller time is re-selected as Δt1 . Δt1 is then put into Eq.(6) as the new iterative time for the k–1 step. The cell i temperature is then k −1 computed to change from Ti ,min to Ti1 . In this second stage the times Δt 2′ and Δt 2′′ when the cell i temperature changes from Ti1 to Tm are calculated using Eqs.(8) and (11), respectively. The larger value is Δt 2 . Δt 2 is put into Eq.(6) as the iterative time for the k–1+ Δt1 step. The cell i temperature is then computed to change from Ti1 to Ti 2 . If

Ti 2 < Tf + ΔT − δ , the smaller time is selected as the

265

new iterative time and the above iterative process is repeated. A similar process is repeated until Tf + ΔT − δ ≤ Ti n < Tf + ΔT . At that moment, n

Δti = ¦ Δtm . m =1

3) Searching for Δtmin

Δtmin represents the minimum time needed for all cells whose temperature reaches the phase change point at the k–1 step. When i=j, the time for the cell i k −1 temperature to go from Ti ,min to Tm is Δti , which

is selected as Δtmin . When ij, the time for the cell j k −1 temperature to go from T j to Tm is Δt j . The

smaller of Δti or Δt j is selected as Δtmin . 4) Selecting Δtmin is calculated for cell i, whose temperature lies in the range between Tf + ΔT − δ and Tf + ΔT , by completing Step 2). Cell i now avoids the “phase change missing” problem. The other cell temperatures lying in the same range at the same time also avoid “phase change missing.” The efficiency is thereby increased. If “phase change overflowing” exists the method is similar. Suppose Tm = Tf − ΔT − δ , cell i, whose temperature lies in the range between Tf − ΔT − δ and Tf − ΔT , avoids “phase change overflowing.” If “phase change missing” and “phase change overflowing” exist together the smaller time is selected after comparing their respective Δtmin values. 5) The iteration continues and Steps 1) to 4) are repeated until the end of the computing time. The calculation of soil thawing is similar to the method mentioned above. Following the above-mentioned method allows the location of the phase change interface to be calculated accurately. “phase change missing and overflowing” is eliminated and the precision of the calculations is improved. This method is an advance over the traditional method.

5

Example analysis

A mathematical model for underground soil cooling has been established. The buried tube is a vertical U-tube heat exchanger. This U-tube is considered equivalent to a differently sized round tube[12] with an equivalent diameter of Deq = 2 D , where D is outer

diameter. For the assumed conditions of the model refer to Reference [13]. The heat transfer model of the temperature around the tube can be treated as axially symmetric. The physical model is shown in Fig. 1.

266

Mining Science and Technology

% “  U

%RUHKROH (TXLYDOHQW URXQGWXEH

δ =0.0005 °C and the fluid Reynolds number is Re=2400. The inlet fluid temperature of the tube during the cold-storing period is shown in Table 1.

)DUILHOG ERXQGDU\

6RLO

 '

[

Table 1

%RWWRP & ERXQGDU\

Inlet fluid temperature of the tube at different times

Inlet fluid temperature (°C)

Fig. 1 Physical model of a vertical U-tube heat exchanger

The two-dimensional unsteady state differential equation for heat conduction in cylindrical coordinates is known (Eq.(5)). In the calculation regions shown in Fig. 1 the AB side of the soil and the AD side of the tube wall are the third boundary condition. The radius of thermal influence, BC, and the bottom boundary, CD, are the adiabatic boundary conditions. The pipe outer diameter is 32 mm. The soil is gravel with a ρd =1800 kg/m3 (density of dry soil), φ =30%,

λf =3.05 W/(m·°C), λu =2.18 W/(m·°C), Cf = Table 2

–2

–5

–5/–12

–5

–2

Daily continuous cold-storing time (h)

12

12

12/12

12

12

Number of days (d)

15

15

40

15

15

Using the completely implicit scheme, the fixed time step method (FTM) or the variable time step method (VTM) were used to find model estimates. When “Phase change missing and overflowing” exists the calculation returns to the previous step and the time step is reduced to 0.2 times the original time step. The BCTM results are shown along with the FTM and VTM results in Tables 2 and 3 and in Figs. 2–4. The objects are the soil close to the U-tube.

Comparison of the results calculated using the FTM, VTM and BCTM models: soil frozen at 0 °C

t (h)

T (°C)

t (h)

T (°C)

tmodify (h)

Tmodify (°C)

BCTM Ƹt1 (s)

2.2 2.6

727.00 727.00

–0.1938 –0.1922

726.20 726.20

–0.0210 –0.0199

726.1208 726.1265

–9.33E-05 –4.97E-04

423.100 20.527

3.0

727.00

–0.1778

726.20

–0.0057

726.1769

–2.55E-04

181.550

1.8 3.4

727.00 727.00

–0.1748 –0.1554

726.20 726.40

–0.0013 –0.0346

726.1938 726.2592

–4.73E-04 –1.06E-04

60.917 235.490

3.8 1.4

727.00 727.00

–0.1280 –0.1217

726.40 726.60

–0.0074 –0.0482

726.3646 726.3981

–3.44E-05 –4.05E-04

372.990 120.560

6.492

4.2

727.00

–0.0976

726.60

–0.0251

726.4865

–1.91E-04

313.780

4.478

X (m)

No.2

1994.8 kJ/(m3·°C) and Cu =2785.2 kJ/(m3·°C). The soil initial temperature is taken from the test results reported in Reference [11]. The time step is 3600 s. The calculation space is divided into 106×17 cells. During calculation Tf =–0.25 °C, ΔT =0.25 °C,

6RLOVXUIDFH $

Vol.19

FTM

VTM

Ƹt2 (s)

Ƹt3 (s)

8.397

3.208

Note: All temperatures shown are higher than 0 °C at 726 hours.

Table 3

Comparison of results calculated using the FTM, VTM and BCTM models: soil thawed at –0.5 °C

t (h)

T (°C)

t (h)

T (°C)

tmodify (h)

Tmodify (°C)

BCTM Ƹt1 (s)

Ƹt2 (s)

Ƹt3 (s)

21.5 21.0 20.6

1693.00 1693.00 1693.00

0.0670 –0.0524 –0.0770

1692.20 1692.20 1692.20

–0.3501 –0.4018 –0.4840

1692.0331 1692.0829 1692.1189

–0.4999 –0.4996 –0.4996

96.192 103.440 107.306

21.735 72.861 18.559

1.223 3.069 3.824

1.0 20.2

1693.00 1693.00

0.2363 –0.0857

1692.20 1692.40

–0.4238 –0.3149

1692.1279 1692.1539

–0.4997 –0.4997

26.862 93.314

4.030

1.451

19.8 19.4

1693.00 1693.00

–0.0716 –0.2275

1692.40 1692.40

–0.4055 –0.3950

1692.1832 1692.2224

–0.4998 –0.4999

104.650 138.938

0.918 2.242

19.0

1693.00

–0.2506

1692.60

–0.4448

1692.2697

–0.4999

166.440

3.750

X (m)

FTM

VTM

Note: All temperatures shown are lower than –0.5 °C at 1692 hours.

In Tables 2 and 3 “phase change missing” (shown in bold-face) of different degrees occurs during the calculations with the FTM and VTM approaches. The FTM’s “phase change missing” is lager than the VTM’s. The FTM has cells 2.2 m deep missing 39% of the phase change latent heat. The cells 1.0 m and 21.5 m deep, whose temperature exceeds 0 °C, miss all latent heat during the thawing process. This error is larger.

The VTM still shows low precision and the cell 20.2 m deep misses 37% of the phase change latent heat during thawing. In general, the BCTM needed to re-iterate the temperature once, and no case did so more than 3 times. All the re-iterated temperatures are within the error range (0 ~ 0.0005 °C). Consequently, the BCTM not only ensures very high precision, but also increases the computing time by only a little. The non-convergence problem of numerical oscillation

YANG Tao et al

An improvement to the simulation of phase-change heat-transfer …

did not occur. If the precision requirements are reduced the time for re-iteration can be reduced. The soil around the tube freezes during the cold storage process. When the system stops storing cold, or the cold-storing temperature rises during the day, parts of the frozen soil thaw. Consequently, a freezing and thawing cycle occurs in the soil. When the FTM or the VTM have “Phase change missing” problems during freezing the calculated temperature is lower and the velocity of the freezing front is faster. When the FTM or the VTM has “phase change overflowing” problems during soil thawing the calculated temperature is higher and the velocity of the freezing front is slower. Because the soil temperature in the non-phase-change region changes much faster than in the phase change region the missing phase change latent heat is greater than the overflowing. This missing latent heat accumulates continuously as time passes leading to a rise in the velocity of the freezing front. The velocities of the freezing front determined in different ways are plotted in Fig. 2. As can be seen,  

)70 970 %&70

     

Fig. 2







  [ P







Time of different depths soil fifth frozen at 0 °C

(a) At t=840 h

Fig. 3

6

267

the FTM and the VTM give a velocity larger than the BCTM. The velocity from the FTM is the fastest and its error is the biggest. These simulations prove the above analysis is right. A laboratory experiment of seasonal soil cold storage was first established in the energy-saving laboratory of Harbin Institute of Technology in October, 2004. The experimental research was continued for a whole year[13–14]. The measured temperatures at depths of 5 m, 10 m, 30 m and 50 m, and the calculated temperatures at these different depths, are seen in Fig. 3. In Fig. 3a the predicted temperatures for the BCTM are smooth and, except at the depth of 50 m, the calculated results are consistent with the test results. The temperature curves from the FTM have larger fluctuations. Temperatures change from 0.9 °C to –0.26 °C rapidly at 12 m depth, which is inconsistent with reality. Additionally, although the time step is reduced, the VTM shows phase change errors. The temperature variation from VTM simulation has irregular fluctuations at the 15 m to 17 m depths, which is also inconsistent with reality. In Fig. 3b the BCTM results are shown to be consistent with the test results. The temperatures at 5 m and 30 m depths calculated by the FTM or the VTM have large differences compared to test results. Temperature curves for FTM or VTM models appear irregular. In these two figures the biggest error between the BCTM and the FTM, or between the BCTM and the FTM, occurs near the freezing front. This is consistent with reality since most of the phase change associated with freezing and thawing occurs near the freezing front. This is where the phase change latent heat is missed by the FTM and the VTM.

(b) At t=1200 h

Temperature curves for soil close to the wall of the well at 840 and 1200 h

Conclusions

In this paper a sensitive heat capacity model is introduced and a new numerical method, which is called the backwards calculating time step method, is presented for solving the phase change heat transfer problem. The two-dimensional problem of the moving boundary is simulated and compared to test results. The temperature fields calculated by the new method are more consistent with measured values.

The method uses an iterative scheme and adjusts the time step automatically. This allows the re-iterated temperature of a cell to fall on the phase change point and so gives more precise rates than traditional methods. It also avoids the non-convergence problem of numerical oscillation during calculations. The method not only can be used for freezing and thawing problem of soil but also may be applied to other complex multi-dimensional problems involving phase change heat transfer.

268

Mining Science and Technology

References [1]

[2]

[3]

[4] [5]

[6]

[7]

Zhou Y, Zhou G Q, Zhang Q. An expanded quasi-static model of ice lens growth during frost heave. Journal of China University of Mining & Technology, 2008, 37(3): 333–338. (In Chinese) He M C, Zhang Y, Guo D M, Qian Z Z. Numerical analysis of doublet wells for cold energy storage on heat damage treatment in deep mines. Journal of China University of Mining & Technology, 2006, 16(3): 278–282. Yu X, Nelson D J, Vick B. Phase change with multiple fronts in cylindrical systems using the boundary element method. Engineering Analysis with Boundary Elements, 1995(16): 161–170. Voller V, Cross M. Accurate solutions of moving boundary problems using the enthalpy method. International Journal Heat Mass Transfer, 1981(24): 545–556. Tanyhienn C, Cheunn F B, Shiah S W. Behavior of the two-phase mushy zone during freeze coating on a continuous moving plate. Journal of Heal Transfer, 2002, 124(1): 111–119. Kong X Q. Time finite difference scheme and its variable step size calculations in solution of non-steady conduction problems by finite element method. Journal of Engineering Thermo Physics, 1982, 3(8): 263–269. (In Chinese) Shang X Y, Zhou G Q, Bie X Y. Numerical simulation improvement of freezing soil’s temperature field. Jour-

[8] [9]

[10] [11] [12]

[13]

[14]

Vol.19

No.2

nal of China University of Mining & Technology, 2005, 34(2): 179–184. (In Chinese) Liu G. Numerical Simulation of Temperature Fields. Chongqing: Chongqing University Press, 1990. (In Chinese) Ye H, He H F, Ge X S, Xu B. The comparative numerical investigations on the melting process of form-stable phase change material using enthalpy formulation method and effective heat capacity formulation method. Acta Energiae Solaris Sinica, 2004, 25(4): 488–491. (In Chinese) Lu J F. Numerical Solution of Partial Differential Equations. Beijing: Tsinghua University Press, 1987. (In Chinese) Tao W Q. Numerical Heat Transfer: Second Edition. Xi’an: Xi’an Jiaotong University Press, 2005. (In Chinese) Gu Y, O’Neal Dennis L. Development of an equivalent diameter expression for vertical U-tube used in groundcoupled heat pumps. ASHARE Trans, 1998, 104(2): 347– 335. Yang T, Zheng M Y. Experimental research and simulation by soil cold storage using natural energy. In: Proceedings of the 22nd International Congress of Refrigeration. Beijing: Golden Info-bridge Publishers, 2007. Liu W. Experimental Research on Natural Energy Air-conditioning System by Soil Cold Storage in Severe Cold Climate [Master dissertation]. Harbin: Harbin Institute of Technology, 2005. (In Chinese)