Study on the influence of water flow on temperature around freeze pipes and its distribution optimization during artificial ground freezing

Applied Thermal Engineering 135 (2018) 435–445

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Study on the influence of water flow on temperature around freeze pipes and its distribution optimization during artificial ground freezing

T

⁎

Shibing Huanga,b, , Yunlin Guoa, Yanzhang Liua,b, Lihua Kea, Guofeng Liuc, Cheng chena a

School of Resources and Environmental Engineering, Wuhan University of Science and Technology, Wuhan, Hubei 430081, China Hubei Key Laboratory for Efficient Utilization and Agglomeration of Metallurgic Mineral Resources, Wuhan 430081, China c School of Highway, Chang’ an University, Xi’an, Shaanxi 710064, China b

H I G H L I G H T S coupled hydro-thermal model was applied for artificial ground freezing. • AThedeveloped of water flow on freezing temperature around freeze pipes was investigated. • A newinfluence combining optimization method of freeze pipes arrangement was proposed. • Three common distribution functions were suitable to determine the optimum positions of freeze pipes. •

A R T I C L E I N F O

A B S T R A C T

Keywords: Artificial ground freezing Coupled hydro-thermal model considering water/ice phase transition Seepage flow Nelder–Mead simplex method Distribution function

Artificial ground freezing is usually used to improve ground and provide temporary support. However, the freezing process is dramatically influenced by water flow. Evidently, water flow with high velocity brings large heat energy that prevents the freezing of porous media around freeze pipes. In this paper, for safety and energysaving, the influence of water flow on freezing process is simulated by a developed coupled hydro-thermal model considering water/ice phase transition and the positions of freeze pipes around circular tunnel are optimized through combining this model with Nelder–Mead simplex method based on COMSOL multiphysics platform. According to the evolution law of freezing band under well-distribution of freeze pipes, three kinds of potential frequently-used probability distribution functions are adopted to reduce control parameters and improve the efficiency of optimization program including normal distribution, Poisson distribution and chi-square distribution. The results show that the proposed combining simulation method is suitable for optimization of freeze pipes arrangement. The freezing time is significantly reduced no matter which one of the above distribution functions is employed through careful design.

1. Introduction Freezing of porous media has been widely studied, either natural or artificial ground freezing [1]. Artificial ground freezing method, providing temporary support and waterproof layer in geotechnical engineering construction, are extensively used in soft rock and soils [2,3]. The fundamental principle in ground freezing is to inject cold energy and convert pore water into ice [4]. When the underwater tunnels pass through the rivers or straits, the process of construction in soft rock or soils usually is very difficult because of inrush disaster. Moreover, water may flow around the underwater tunnels before excavation because of the difference of water head. Artificial ground freezing technique is very suitable to improve the soil strength and prevent water flow into

⁎

the excavated area in this case. To meet the demand of frozen soil strength, the thickness of frozen arch should be greater than 1.5 m [5]. However, water flow brings much thermal energy against the formation of freezing arch and changes the freezing path. The required freezing time of forming a specified frozen arch is also considerably influenced by high seepage flow. It has been investigated the area between two freeze pipes may not freeze if the flow velocity of water exceeds 1–2 m/ d in high permeability soils or fractured media [6]. Moreover, when the temperature drops below freezing point, partial water in porous media gradually freezes but some liquid water is also existed [7,8]. Therefore, the freezing process of porous media is very complicated but important to conduct artificial freezing construction, which is associated with the coupled hydro-thermal action under low temperature.

Corresponding author at: School of Resources and Environmental Engineering, Wuhan University of Science and Technology, Wuhan, Hubei 430081, China. E-mail address: [email protected] (S. Huang).

https://doi.org/10.1016/j.applthermaleng.2018.02.090 Received 23 September 2017; Received in revised form 10 February 2018; Accepted 25 February 2018 Available online 27 February 2018 1359-4311/ © 2018 Elsevier Ltd. All rights reserved.

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

Many freezing models of coupled heat and water flow, freezing experiments and related numerical studies on porous media have been reported in the past decades [9–12]. However, the influence of high seepage flow on freezing process is rarely considered [13]. Besides, at present, there are a lot of researches on freezing construction technology and case analysis about artificial ground freezing, in which freeze pipes are arranged evenly or empirically [14–16]. Nevertheless, the freezing rates between two adjacent freeze pipes are unequal in the presence of seepage flow and freeze pipes may satisfy a certain distribution law which needs to be further studied [17]. To our knowledge only Marwan et al. [4] has tried to optimize the positions of freezing pipes by Ant Colony Optimization and it results in a significant decrease of freezing time. Due to the special freezing process of porous media under high seepage velocity, it is necessary to study the freezing process and propose corresponding optimization method during artificial ground freezing for safety, energy saving and high efficient construction. In this paper, the Nelder–Mead simplex method is adopted to optimize the positions of freeze pipes through combining it with a developed coupled hydro-thermal model in order to derive a required frozen arch as fast as possible. The Nelder-Mead simplex method is an efficient and derivative free optimization algorithm to find the minimum or maximum value of an objective function, which has been widely applied in many fields, including engineering materials, geotechnical and hydrologic engineering [18,19]. It is relatively simple and suitable for not too many control parameters, in which only a numerical evaluation of the objective function is required [20]. The basic idea for the simplex algorithm from geometry is shown in Fig. 1. In the three dimensional space, a simplex is a special tetrahedron determined by four points (F1, F2, F3 and F4) and their connected lines. The objective function is estimated at every point. The highest point, where the objective function is largest (e.g. Point F3), will be perpendicularly mirrored against the opposite plain segment to a lower point, which is also accompanied with an expansion or contraction to modify step size in order to reach the optimization valley floor. A termination criteria should be given to stop this optimization procedure, generally including the maximum number of reflections or a tolerance for critical variables. Similarly, it can be extended to the N dimension in which the simplex is a special polytope of N + 1 vertices. Here a coupled hydro-thermal model considering the influence of water/ice phase transition, which is firstly studied by Tan et al. [22], is proposed and validated by a previous large-scale experiment of artificial ground freezing considering the effect of water flow on the freezing process in Section 2. Then an optimization method of freeze pipes arrangement, combining the coupled hydro-thermal model with an

efficient optimization algorithm (Nelder-Mead simplex method), is proposed and applied in a circular tunnel constructed by artificial ground freezing method in Section 3. The freezing temperature around freeze pipes and the minimum freezing time under different seepage flow conditions is investigated before and after optimization in Section 4. In Section 5, some issues about the topic are discussed. Finally, some significant conclusions are drawn in the last section. 2. Coupled hydro-thermal model considering water/ice phase transition Thermal transfer and water flow during freezing around freeze pipes should be considered when investigating the influence of water flow on freezing temperature around freeze pipes during artificial ground freezing. Therefore, two crucial functions should be presented, heat conduction equation with phase change deduced from energy conservation and continuity equation considering water flow deduced from water mass conservation. 2.1. Governing equations 2.1.1. Basic assumptions Considering the actual physical process, some basic assumptions are introduced as below: (1) The rocks or soils are saturated, homogeneous and isotropic porous media; (2) the evaporation process of water is ignored, and the Darcy’s law is suitable to describe the groundwater flow in porous media; (3) the heat conduction in freezing porous media satisfies Fourier’s law. 2.1.2. Heat conduction equation A fully coupled hydro-thermal model has been firstly proposed by Tan et al. [22]. The heat conduction equation deduced from energy conservation in freezing porous media is expressed as

Cv

∂T + ρl cl vl ·(∇T ) + ∇ ·(−λ e ∇T ) = 0 ∂t

(1)

where ρ and c are density and specific heat capacity, respectively. Subscript s, l and i represent solid matrix, water and ice, respectively. T is temperature. vl is the seepage velocity vector of water. Water density is the function of pressure and temperature:

ρl = ρl0 [1 + αT (T −T0) + βp (pl −p0 )]

(2)

where ρl0 is the density corresponding to the initial pressure p0 and initial temperature T0 ; αT is the thermal expansion coefficient of water related with temperature; βp is water compressibility; pl is water presp0 = 101 kPa T0 = 20 °C, and there are sure. When αT = (−9T −80) × 10−6 /°C and βp = 5 × 10−6 /Pa. The change of water density with temperature when p0 = 101 kPa is shown in Fig. 2. λ e is the effective thermal conductivity of porous media, which depends on the thermal conductivity of its components (solid, unfrozen 1000 999

ρl (kg/m3)

998 997 996 995 994 993 -25

-20

-15

-10

-5

0

5

10

15

20

25

T (°C) Fig. 1. Nelder-Mead simple for three optimization parameters [21].

Fig. 2. The change curve of water density with temperature [23].

436

30

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

water and ice), their volumetric fractions and the spatial distribution of its components. The exponential weighted mean model is usually used to describe its change rule [24,25]:

λ e = λs1 − n λln·wu λin·(1 − wu)

(3)

where n is porosity and λ is thermal conductivity.wu is the unfrozen water content, which could be expressed as an exponential function of temperature [26]:

wu = e M (T − Tm)

(4)

where M is a parameter related with the distribution of pore radius. Tm = 0 °C is the freezing point of bulk water.Cv is the equivalent volumetric thermal capacity, which can be expressed as follows:

Cv = (1−n) ρs cs + n (1−wu ) ρi ci + nwu ρl cl + nρl ℓ·

∂wu ∂T

(5)

where ℓ is the latent heat of water/ice phase transition; c is specific heat.

Fig. 3. Plan view at measuring plane with location of thermistor.

located in the middle line. Several thermistors of type Pt 100 are also placed along two middle lines (ML1 and ML2) inside the porous media to monitor the freezing temperature field, which are indicated by black solid points (the detailed coordinate positions can be found in the original literature). This large-scale physical simulation model can be simplified to a plane problem as is shown in Fig. 3. The unfrozen water content of this porous material under different freezing temperature was measured using Time-Domain Reflectometry (TDR) in the experiment. Then, we can derive the unknown parameters M = 3.93 from Eq. (4) by fitting the experimental results as is shown in Fig. 4. Obviously, unfrozen water content can be well expressed by the exponential function of temeprature. It should be noted that the calculated freezing point of 0.1 °C is slightly greater than theoretical value 0 °C, which may be caused by measurement error. The other physical and thermodynamics parameters for porous material and water/ice media are listed in Tables 1 and 2, respectively. The segregation potential of this medium sand is very small, which may be negligible [33]. Two different seepage boundary conditions (Initial seepage velocities: v0 = 0 m/d and v0 = 2 m/d) were investigated in this experiment. Thermal boundary conditions of freeze pipes and water inlet gauged by the temperature sensors are given in Fig. 5. Numerical analysis is carried out using the developed coupled hydro-thermal model under the same conditions. The calculated results of temperature at all measured positions are in good agreement with the experimental values after different freezing time as is shown in Fig. 6. Obviously, the freezing rates in porous media without water seepage are much faster. After freezing of 45 h, all the measured points on Line ML2 are frozen with the absence of water flow while many measured points on Line ML2 but

2.1.3. Continuity equation The continuity equation of water/ice has been derived by Tan et al. [22] as below:

∂ [ρi (1−wu ) n] ∂ (ρl wu n) + ∇ ·(ρl vl) = 0 + ∂t ∂t

(6)

Assume water seepage in porous media satisfies Darcy’s law. Considering the effect of segregation potential, water migration in freezing porous media is

vl = −

kr k (∇pl −ρl g)−SP0 ∇T μl

(7)

ks 0 ⎤ is the intrinsic permeability matrix, in which ks is where k = ⎡ ⎢ ⎦ ⎣ 0 ks ⎥ 0 the permeability of saturated porous media before freezing; g = ⎡ g ⎤ is ⎣ ⎦ the gravity acceleration vector.SP0 is segregation potential, which is a function of average suction in the frozen fringe, the applied pressure on the porous media and the rate of cooling [27–29]. For example, the average value of segregation potential of Calgary silt is approximately 2.3 × 10−9 m2/°C/s in segregation interval −0.3 to −0.1 °C [30]. μl is the viscosity of liquid water which could be expressed as a function of freezing temperature [23]: 1808.5 ⎞ μl = 2.1 × 10−6exp ⎛ ⎝ 273.15 + T ⎠

(8)

kr is the relative permeability which could be written as a function of water saturation [31]: kr =

Sl [1−(1−Sl1/ κ )κ ]2

(9)

(

)

e M (T − Tm) [1− 1−e M (T − Tm)/ κ κ ]2

Unfrozen water content

kr =

Experiment values Fitting curve

1.0

where κ is material constant and Sl is water saturation. The water saturation is equal to unfrozen water content during freezing, then we can derive the final expression of relative permeability as below: (10)

2.2. Model validation Pimentel et al. [32] had conducted a large-scale laboratory test on artificial ground freezing under different initial seepage flow velocities from 0 m/d to 2.1 m/d. Medium sand was filled in a watertight box with inner dimensions 1.3 m × 1.0 m × 1.2 m to simulate porous material in this experiment. The box has a stable water inlet at left and water outlet at right. All the surfaces have thermal insulation with at least 18 cm of foam to create adiabatic boundary conditions. Three freeze pipes with diameter of 0.041 m perpendicular to the seepage-flow direction are

0.8 0.6

wu = e3.93(T −0.1)

0.4 0.2 0.0 -2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

Temperature (°C) Fig. 4. Measured and calculated unfrozen water content against temperature.

437

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

same as that in the above experiment as is shown in Table 1.

Table 1 The values of thermo-mechanical parameters of experiment material [32]. n

ρs (kg/ m3)

λs (W/m/ °C)

cs (J/kg/ °C)

ks (m/s)

SP0 (m2/ °C/s)

M (/°C)

κ

0.41

1571

4.3

816

0.0001

0

3.93

0.5

3.2. The influence of water seepage on freezing temperature The development of frozen arch (T ⩽ −3 °C) under well-distributed pipes before optimization could be simulated by solving the developed coupled hydro-thermal model under related boundary conditions above using finite element method as is shown in Fig. 8. The initial flow velocity of water v0 is 1.5 m/d in this case. It shows that a long frozen band is formed to connect the freeze pipes arranged along the flow direction on the top of designed arch after freezing of 10 days. However, there is no freezing intersecting band produced on the upstream or downstream side. This frozen band grows more quickly on the downstream side along with freezing time because of the resistance of heat inflow by freeze pipes on the upstream side as is shown in Fig. 8(b). Water flows much faster through the unfrozen area of adjacent freeze pipes on the upstream and downstream sides due to the freezing of the top side before a complete frozen arch is formed (larger arrows represent higher velocity in Fig. 8). After freezing of 30 days, a complete frozen arch arises. Then, the seepage flow of water insides excavation zone stops because of this waterproof arch. Water keeps flowing to the downstream sides through the unfrozen area above the frozen arch. However, the frozen arch does not occur in the designed ring and shifts to the left. The width of frozen band on the upstream side is much smaller than that on the downstream owing to the influence of water flow. Therefore, it is necessary to optimize the positions of freeze pipes in order to generate a designed frozen arch using significantly less freezing time.

Table 2 The values of calculation parameters of water/ice.

l (J/kg)

ρi (kg/ m3)

cl (J/kg/ °C)

ci (J/kg/ °C)

λi (W/m/ °C)

λl (W/m/ °C)

334.88 × 103

917

4200

2100

2.2

0.6

far from freeze pipes are unfrozen under water flow of 2 m/d due to the inflow of heat. Therefore, the seepage flow will hinder freezing around freeze pipes and result in the increase of freezing time when a frozen band connecting freeze pipes is created. Actually, the calculation accuracy of the proposed coupled hydro-thermal model under low temperature are greatly influenced by the critical coupling parameters, mainly including the unfrozen water content and the effective thermal conductivity of freezing porous media. Through the comparisons between the calculated and measured temperature as is shown in Fig. 6, those two important parameters can be well expressed by the exponential weighted mean model and the exponential function, respectively. 3. Optimization of artificial ground freezing considering seepage flow

3.3. Optimization of freezing pipes arrangement The positions of N freeze pipes could be determined by the polar coordinates (ri,θi ) . 2 × N unknown parameters need to be optimized if the distribution law of freeze pipes is not considered. Therefore, some potential distribution functions may be adopted to reduce the numbers of optimized parameters and improve the efficiency of optimization. The radial coordinates of the freeze points (the centers of those freeze pipes) could be written as (Fig. 9)

3.1. Model description Consider artificial ground freezing for circular tunnel excavation under seepage flow similar to the case studied by Marwan et al. [4], which is originally performed by Ziegler et al. [17], as is shown in Fig. 7. Eighteen brine freeze pipes with steady cooling temperature of −30 °C and diameter of 0.1 m are arranged to freeze soils around designed section of this tunnel in order to form a closed frozen arch with expected thickness of 1.5 m. The frozen arch is cooled to −3°C for ensuring the strength and stability [33,35]. A stable water flows from the left to right insides saturated porous media with initial temperature of 10 °C. According to the scale and characteristic of freezing project, half of the study area is taken to build a numerical model (width: 30 m, depth: 10 m). The related parameters of porous media are chosen as the

(11)

ri = r0 + δr (θi )

where r0 is the radius of center circle of designed arch; δr (θi ) is the radially shift distance of freeze points after optimization. δr (θi ) may be assumed to be linear variation with θi on the downstream side and constant on the upstream side, which could be expressed as below [4]:

25

25

20

20

Th (Freeze pipe)

15

15

Th (Freeze pipe)

5

5

Tw0 (Water inlet)

T (°C)

10

T (°C)

10

0 -5

0 -5

-10

-10

-15

-15

-20

-20

-25

0

5

10

15

20

25

30

35

40

-25

45

0

5

10

15

20

25

t (h)

t (h)

(a) v0 = 0 m/d

(b) v0 = 2 m/d

Fig. 5. Thermal boundary conditions of freeze pipes and water inlet.

438

30

35

40

45

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

30

Experiment results

25

T10 T22 T10 T22

Developed model

20 15

T11 T36 T11 T36

20

T13 T37 T13 T37

10

0

T30 T35 T30 T35

-5

-5

-10

-10

-15

-15

-20

-20

-25 0

10

20

30

40

t (h)

50

60

70

80

90

0

10

20

30

30

Experiment results

25

Developed model

20

40

50

60

70

80

90

100

t (h)

(b) v0 =0 m/d (ML2)

(a) v0 =0 m/d (ML1)

T10 T22 T10 T22

T11 T36 T11 T36

25

T13 T37 T13 T37

Experiment results T31 Developed model T31

20 15

T (°C)

15

T (°C)

T29 T34 T29 T34

0

T (°C)

T (°C)

5

10

T28 T33 T28 T33

T29 T34 T29 T34

T30 T35 T30 T35

10 5

5 0

0

-5

-5

-10

T28 T33 T28 T33

5

10

-25

Experiment results T31 Developed model T31

15

0

5

10

15

20

25

30

35

40

45

-10

50

0

5

10

15

20

t (h)

25

30

35

40

45

50

t (h)

(c) v0 = 2 m/d (ML1)

(d) v0 = 2 m/d (ML2)

Fig. 6. Calculated and measured temperature against freezing time.

The other position parameter need to be characterized is the angular coordinate θi . There is

θi =

θi ⩽ π /2 ⎧ d2 2θ ⎨ (d1−d2) π i + 2d2−d1 θi > π /2 ⎩

0 i=1 θi − 1 + δθi − 1 I > 1

(13)

where δθi − 1 is the difference of angular coordinate between the i th and (i − 1) th freeze points. The positions of freeze pipes after optimization in polar coordinate system are shown in Fig. 9 (red1 solid points).As we know, a complete frozen band between two adjacent freeze pipes must occur on the top of designed arch firstly and extend to both sides under well-distribution freeze pipes as is shown in Fig. 8. Besides, the freezing of porous media on the upstream side is much harder due to water flow. Therefore, the following three condition should be satisfied when rearranging the freeze pipes: (1) Freeze pipes should be more compactly rearranged on the upstream side to reduce freezing time of total designed arch; (2) δθi − 1 between two adjacent freeze pipes should increase with the increasing of θ from left to right until reaches the maximum value; (3) the maximum δθi − 1 should be located on the position after the fourth freeze pipe where i ≥ 5. By the above analysis, three potential and common distribution functions are adopted to rearrange the freeze pipes, including normal distribution, Poisson distribution and chisquare distribution, because it’s hard to determine which one is much better before optimization.

Fig. 7. Calculated model of artificial ground freezing.

δr (θi ) =

{

(12)

where d1 and d2 are the maximum radially shift distance of freeze points on the downstream side and on the upstream side, respectively. Considering the designed thickness of frozen arch (1.5 m) and the radial development of frozen band around freeze pipes as is shown in Fig. 8, they can be set as −0.5 ⩽ d1 ⩽ 0 and 0 ⩽ d2 ⩽ 0.5, respectively. It means that all the freeze pipes should be move from downstream toward upstream side as is shown in Fig. 9.

1 For interpretation of color in ‘Fig. 9’, the reader is referred to the web version of this article.

439

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

(a) 10 d

(b) 20 d

(c) 30 d

Fig. 8. Frozen arch by well-distributed pipes (v0 = 1.5 m/d, including the velocity vector of water flow).

(1) Normal distribution The normal distribution is a very common continuous probability distribution. The probability density function is

(x −μ)2 ⎞ 1 exp ⎛− 2σ 2 ⎠ 2π σ ⎝

f (x ,μ,σ ) =

⎜

⎟

(14)

where μ and σ are distributed parameters representing mathematical expectation and standard deviation, respectively. As δθi − 1 should gradually increase from upstream side to the maximum value on the downstream side. It may be expressed as Fig. 9. Positions of freeze pipes after optimization in polar coordinate system.

δθi − 1 = π

f (i−1,μ,σ ) 10

∑i = 2 f (i−1,μ,σ )

(15)

Obviously, when 〈i−1−μ〉 = 0 , δθi − 1 gets the maximum values. Therefore, the maximum δθi − 1 occurs between the 〈μ〉th and (1 + 〈μ〉)th freeze pipe. Where 〈〉 is the rounding symbol. It has been validated that the maximum δθi − 1 should occur after the fourth freeze pipe where 5 ⩽ i ⩽ 10 , because the frozen band starts from the fourth freeze pipe as is shown in Fig. 8. Therefore, we have 4 ⩽ μ ⩽ 9 and 4 ⩽ σ ⩽ 9. (2) Poisson distribution Poisson distribution is an important discrete probability distribution, which can express the probability of a given number of events occurring in a fixed interval space if these events happen with a known average rate. The probability density function can be expressed as

g (x ,λ ) =

(16)

where x > 0 ; λ represents the mathematical expectation and standard deviation. According to the properties of probability density function of Poisson distribution, when x = λ , g (x ,λ ) gets the maximum value. As we know the maximum δθi − 1 should be located between the 〈4 + k〉th and 〈5 + k〉th freeze pipes where k ⩾ 0 . Then we can define δθi − 1 as below

Fig. 10. Flowchart of the optimization.

(a) 5 d

λ x e−λ x!

(b) 15 d Fig. 11. Frozen arch after optimization adopting normal distribution (v0 = 1.5 m/d).

440

(c) 25 d

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

(a) 5 d

(b) 15 d

(c) 24 d

Fig. 12. Frozen arch after optimization adopting Poisson distribution (v0 = 1.5 m/d).

(a) 5 d

(b) 15 d

(c) 25 d

Fig. 13. Frozen arch after optimization adopting chi-square distribution (v0 = 1.5 m/d).

Completed freezing time tf (d)

120

80

Normal distribution Poisson distribution Chi-square distribution Uniform distribution Fitting curve before optimization Fitting curve after optimization

60

t f = 16 + 2.53e

100

distribution of a sum of the squares of q independent standard normal random variables, which is widely used in inferential statistics. The probability density function is

p (x ,q) =

1.79 v0

δθi − 1 = 1.0

1.5

2.0

2.5

3.0

Initial seepage velocity v0 (m/d)

g (λ−(5 + k ) + i,λ ) 10

∑i = 2 g (λ−(5 + k ) + i,λ )

(17)

Obviously, when λ−〈5 + k〉 + i = λ , δθi − 1 reaches the maximum value. Then we have i = 5 + 〈k〉. k is a unknown control parameter need to be optimized. Then we can derive 0 ⩽ k ⩽ 5, because there is 5 ⩽ i ⩽ 10 when δθi − 1 takes the peak value. Another unknown parameter λ is from 10 to 70 in order to ensure that the freeze pipes can be sufficiently intensive layout on the upstream. For convenience, λ is replaced by 10 λ to reduce the optimization interval. Then, we have

δθi − 1 =

g (10λ−(5 + k ) + i,10λ ) 10

∑i = 2 g (10λ−(5 + k ) + i,10λ )

p (q−(7 + s ) + i,q) 10

∑i = 2 p (q−(7 + s ) + i,q)

(20)

From Eq. (20), δθi − 1 reaches the maximum values when q−〈7 + s〉 + i = q−2 . Then we can have i = 5 + 〈s〉. s and q are unknown control parameters need to be optimized. Similarly, there is 0 ⩽ s ⩽ 5 for 5 ⩽ i ⩽ 10 at the peak value of δθi − 1. According to characteristics of probability density function, set 10 ⩽ q ⩽ 24 to ensure that δθi − 1 could be smaller enough on the upstream sides. Therefore, the ranges of optimized control parameters are 0 ⩽ s ⩽ 5 and 10 ⩽ q ⩽ 24 . From the above analysis, no matter which distribution function is chosen to determine the optimal positions of freeze pipes, there are four control parameters need to be optimized for rearranging the freeze pipes by Nelder-Mead simple method. The basic idea of determining the optimal control parameters by the Nelder-Mead simple method including the following steps: (1) Several initial values of control parameters are inputted and the corresponding values of objective function can be calculated; (2) the decline direction of objective function is determined by evaluating the values of objective function; (3) new control parameters are chosen on this decline direction and the corresponding new values of the objective function are derived; (4) by comparing the new with the original values of the objective function, a new decline direction can be determined. Therefore, during any optimization process, new values of the control parameters are chosen and optimized by Nelder-Mead simple method. New parameter values will

Fig. 14. Comparison of completed freezing times versus initial seepage velocity.

δθi − 1 =

(19)

q

t f = 16 + 0.42e 2 v0

0.5

()

x

()

20

0.0

q

x 2 − 1e− 2

where x > 0 ; Γ 2 denotes the Gamma function. Similarly, when x = q−2 , p (x ,q) takes the maximum value because x is equal to the mode. As the same, the maximum δθi − 1 is between the 〈4 + s〉th and 〈5 + s〉th freeze pipes where s ⩾ 0 . Then we can define δθi − 1 as below:

40

0

1 q q 22Γ 2

(18)

where 1 ⩽ λ ⩽ 7 and 0 ⩽ k ⩽ 5. (3) Chi-square distribution The chi-squared distribution with q freedom degrees is the 441

Applied Thermal Engineering 135 (2018) 435–445

Optimal solution set

d1

d2

Error solution set

d1

d2

0.4

μ

0.6

9 8 7 6 5 4 3

σ

S. Huang et al.

9 8 7 6 5 4 3

d1(d2) (m)

0.2 0.0 -0.2 -0.41

-0.4

-0.50

-0.6 0

10

20

30

40

n

50

60

70

80

90

0

10

20

30

d1

d2

5

Error solution set

d1

d2

4

60

0

10

20

30

40

n

70

80

σ σ

Optimal solution set Error solution set

50

60

70

90

80

Optimal solution set Error solution set

90

k k

k

3 2

0.2

d1(d2) (m)

50

6 Optimal solution set

0.4

1 0

0.0

-0.38 -0.4

λ

-0.2

-0.50

-0.6 0

10

20

30

40

50

n

0.6

60

70

80

7 6 5 4 3 2 1 0

0

0

10

10

20

20

30

30

(b) Poisson distribution

Optimal solution set

d1

d2

Error solution set

d1

d2

s

0.4 0.2

d1(d2) (m)

40 n

(a) Normal distribution 0.6

μ μ

Optimal solution set Error solution set

0.0

40 n

50

40

50

n

q

-0.40

-0.4

-0.50

20

30

40

n

50

60

70

80

60

70

80

0

10

20

30

40

50 n

60

70

s s

80

Optimal solution set Error solution set

90 q q

20

16 10

λ λ

18

-0.6 0

80

Optimal solution set Error solution set

22

-0.2

70

Optimal solution set Error solution set

6 5 4 3 2 1 0

24

60

90

0

10

20

30

40

n

50

60

70

80

90

(c) Chi-square distribution Fig. 15. Determination of the optimum parameters (v0 = 1.5 m/d).

where t is the freezing time. x is the vector of the control parameters. R is the bound set of the control parameters which have been derived above. Tfa is the maximum temperature in the expected frozen arch. Therefore, the less the completed freezing time, the better the freezing quality.

cause rebuilt of geometry model and regeneration of FEM model until the optimal solution is obtained or the maximum number of objective evaluations is reached. This whole optimization procedure is shown in Fig. 10. In the optimization process, the objective function is the freezing time calculated by FEM model when a closed frozen arch (T ≤ −3°C) with expected thickness of 1.5 m is created, which can be expressed as

min t (x) subject to Tfa ≤−3 °C X⊂R

(21) 442

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

method combining with any one of the three distribution functions mentioned above. It means that the Nelder-Mead simplex method is very effective for optimization of freeze pipes arrangement and the chosen distribution functions (normal distribution, Poisson distribution and chi-square distribution) can well characterize the optimal distribution law of freeze pipes under water flow condition by careful design.

Table 3 The most optimum group of parameters for arranging freeze pipes. Normal distribution

v0 (m/d) 0.5 1 1.5 2 2.5

d1 (m) −0.15 −0.19 −0.41 −0.49 −0.49

d2 (m) 0.19 0.22 0.24 0.19 0.18

μ 6.20 6.67 6.82 8.58 7.10

σ 7.59 6.48 5.75 6.51 5.93

Poisson distribution

v0 (m/d) 0.5 1 1.5 2 2.5

d1 (m) −0.23 −0.21 −0.38 −0.43 −0.49

d2 (m) 0.24 0.25 0.25 0.21 0.17

k 2.67 3.10 4.00 4.36 4.53

λ 3.74 3.84 4.00 4.00 4.13

v0 (m/d) 0.5 1 1.5 2 2.5

d1 (m) −0.13 −0.16 −0.40 −0.45 −0.48

d2 (m) 0.20 0.23 0.27 0.28 0.17

s 1.08 1.05 1.67 2.00 0.94

q 21.82 23.52 20.00 20.90 17.06

Chi-square distribution

4.3. The optimum positions of freeze pipes Actually, the minimum freezing time corresponds to several groups of values of control parameters after optimization no matter which distribution function is adopted. Fig. 15 has shown that many combinations of values of control parameters may result in the minimum freezing time of designed arch. All of those values of the control parameters belong to the optimal parameter set, which are represented by the solid points in Fig. 15. The minimum freezing time cannot be obtained by the other parameter set, so those values are attributed to error parameter set, which are represented by the hollow points in Fig. 15. Evidently, the control parameters will converge to be stable values corresponding to the minimum freezing time after dozens of optimization steps. However, in order to reduce the excavation and disturbance area of frozen zone, the excess frozen zone beyond designed frozen arch should be located outside the excavation section as far as possible. Actually, the excess frozen zone on the upstream side is small and there is a little different between the maximum and the minimum d2 of the optimal parameter set as shown in Figs. 11–13 and 15, respectively. Therefore, the minimum d1 of the optimal parameter set may be adopted to reduce the frozen zone in the excavation area on the downstream side. Based on the above principle, the most optimum group of parameters under different flow velocities are derived and listed in Table 3. It is clear that the values of common parameters (d1 and d2) from different distribution functions are very close under the same flow condition. Because the radial movement of freeze pipes is decided by Eq. (12) no matter which distribution function is adopted. Compared with the optimization method raised by Marwan et al. [4], in which the Ant Colony optimization algorithm is applied to search the optimal control parameters, the proposed combining optimization method in this paper have much better stability and reliability, including the control parameters and the minimum freezing time. Moreover, the most optimum d1 gradually increases with the increasing of water flow velocity. The influence of water flow velocity on the most optimum positions of freeze pipes can be seen in Fig. 16. Obviously, the freeze pipes should be increasingly closely arranged on the upstream side and the total freeze pipes gradually move to left as water flows faster. It should be noted that our model is limited to the initial flow velocity of water less than 3 m/d, which is the nature flow velocity of groundwater inside porous media in most cases. When the water flows much fast, there may be turbulence flow, which is beyond the scope of this research. The freezing of porous media under turbulence flow of water need to be researched in the future, including the freezing process of fractured media under much high seepage velocity.

Fig. 16. The most optimum positions of freeze pipes using normal distribution.

4. Results 4.1. Temperature around freeze pipes In this section, the optimization simulation of ground freezing under different initial flow velocities from 0 m/d to 2.5 m/d has been performed. The freezing process of designed arch under seepage flow of 1.5 m/d after optimization are shown in Figs. 11–13. After freezing of 5 days, the maximum frozen zone arises around the sixth or the seventh freeze pipe on the top because of the prevention of heat inflow by the fourth and fifth freeze pipes (the number of freeze pipes increases along counterclockwise). It is encouraging that the frozen bands occur not only on the top but also on the corners of designed arch because of the space reduction of freeze pipes on the corners and extend to the middle until a long and complete frozen arch is formed. Therefore, it is different from the growth process of frozen band with well-distributed pipes. Moreover, the optimal freezing process is very close to each other no matter which distribution function is adopted. As a results, the freezing process is faster and a thicker complete frozen arch is created after about 25 days after optimization.

5. Discussion Because of the injection of heat by water flow, the freezing process of designed arch during artificial ground freezing is much harder, especially the freezing of porous media on the upstream sides. In order to generate a complete frozen arch within limited time, it necessary to improve freezing intensity on the upstream sides or adding straight-row freezing pipes to cut off water flow [36]. Optimization of freeze pipes arrangement is an energy-saving and effortless technique. Therefore, in this paper an optimization method to rearrange freeze pipes by combining finite element analysis of coupled hydro-thermal model with the Nelder–Mead simplex method is proposed. As we known, the release of latent heat caused by water/ice phase transition will prevent the

4.2. The minimum freezing time What’s more, the minimum freezing time of designed arch exponentially increases with the increase of water flow velocity before and after optimization. A desirable prediction function for estimating the minimum freezing time has been derived as is shown in Fig. 14. The completed freezing time of designed frozen arch increases faster under well-distributed pipes, which is about 108 days in the presence of seepage flow of 2.0 m/d. However, only 38 days is needed to create a designed frozen arch under the same seepage flow condition after optimizing the positions of freeze pipes by the Nelder-Mead simplex 443

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

functions is adopted to arrange freeze pipes including normal distribution, Poisson distribution and chi-square distribution. (5) Our proposed method is limited to optimize freeze pipes around circular section at present. More complicated distribution functions need to be established if this optimization method is applied to the other section shapes when using artificial freezing method.

freezing process [37]. Therefore, the influence of latent heat of water/ ice phase transition on freezing is also considered. For the sake of reducing the control parameters and improving calculation efficiency, the optimum position of freeze pipes is assumed to satisfy a particular distribution function including normal distribution, Poisson distribution and chi-square distribution. Then, only four control parameters are chosen to be optimized and related ranges of values are determined by analysis. It should be noted that there may be some other reasonable distribution functions can be used to promote optimization process. We cannot try one by one owing to the limited time and so many probability distribution functions. Therefore, three commonly-used distribution functions are employed in this paper by careful design. As a result, the completed freezing time of designed arch is greatly reduced by adopting this combining simulation method. However, in this paper there are still two issues need to be addressed in the future research:

Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant No. 41702291), the Natural Science Foundation of Hubei Province (No. 2015CFA142), the Hubei Key Laboratory for Efficient Utilization and Agglomeration of Metallurgic Mineral Resources (2017zy005). References

(1) Optimization of freeze pipes arrangement around non-circular cross section

[1] O.B. Andersland, B. Ladanyi, Frozen Ground Engineering, John Wiley & Sons, Inc., 2004. [2] G. Russo, A. Corbo, F. Cavuoto, et al., Artificial ground freezing to excavate a tunnel in sandy soil, measurements and back analysis, Tunn. Undergr. Space Technol. 50 (2015) 226–238. [3] E. Pimentel, S. Papakonstantinou, G. Anagnostou, Numerical interpretation of temperature distributions from three ground freezing applications in urban tunnelling, Tunn. Undergr. Space Technol. 28 (2012) 57–69. [4] A. Marwan, M.M. Zhou, M.Z. Abdelrehim, et al., Optimization of artificial ground freezing in tunneling in the presence of seepage flow, Comput. Geotech. 75 (2016) 112–125. [5] X.D. Hu, S.J. Deng, Ground freezing application of intake installing construction of an underwater tunnel, Proc. Eng. 165 (2016) 633–640. [6] M. Vitel, A. Rouabhi, M. Tijani, et al., Thermo-hydraulic modeling of artificial ground freezing: application to an underground mine in fractured sandstone, Comput. Geotech. 75 (2016) 80–92. [7] E.J.A. Spaans, J.M. Baker, The soil freezing characteristic: its measurement and similarity to the soil moisture characteristic, Soil Sci. Soc. Am. J. 60 (1996) 13–19. [8] Y.M. Lai, W. Pei, M. Zhang, et al., Study on theory model of hydro-thermal–mechanical interaction process in saturated freezing silty soil, Int. J. Heat. Mass. Trans. 78 (2014) 805–819. [9] K. Hansson, J. Šimůnek, M. Mizoguchi, et al., Water flow and heat transport in frozen soil, Vadose Zone J. 3 (2014) 693–704. [10] J.M. McKenzie, C.I. Voss, D.I. Siegel, Groundwater flow with energy transport and water–ice phase change: numerical simulations, benchmarks, and application to freezing in peat bogs, Adv. Water Resour. 30 (2007) 966–983. [11] W.B. Yu, W.B. Liu, Y.M. Lai, et al., Nonlinear analysis of coupled temperatureseepage problem of warm oil pipe in permafrost regions of Northeast China, Appl. Therm. Eng. 70 (2014) 988–995. [12] W. Song, Y. Zhang, B. Li, et al., A lattice Boltzmann model for heat and mass transfer phenomena with phase transformations in unsaturated soil during freezing process, Int. J. Heat Mass Trans. 94 (2016) 29–38. [13] M. Vitel, A. Rouabhi, M. Tijani, et al., Modeling heat and mass transfer during ground freezing subjected to high seepage velocities, Comput. Geotech. 73 (2016) 1–15. [14] X.D. Hu, S.J. Deng, H. Ren, In situ test study on freezing scheme of freeze-sealing pipe roof applied to the gongbei tunnel in the Hong Kong-Zhuhai-Macau bridge, Appl. Sci. 7 (2017) 1–13. [15] Q.X. Yan, Y.J. Xu, W.B. Yang, et al., Nonlinear transient analysis of temperature fields in an AGF project used for a cross-passage tunnel in the Suzhou Metro, KSCE J. Civ. Eng. (2017) 1–11, http://dx.doi.org/10.1007/s12205-017-1118-4. [16] Y.S. Kang, Q.S. Liu, Y. Cheng, et al., Combined freeze-sealing and New Tubular Roof construction methods for seaside urban tunnel in soft ground, Tunn. Undergr. Space Technol. 58 (2016) 1–10. [17] M. Ziegler, C. Baier, B. Aulbach, Numerical simulation of artificial ground freezing applications in tunnelling, in: Computational Methods in Tunneling (EURO: TUN 2009) 2009 Ruhr University Bochum Aedificatio, 2009, pp. 343–350. [18] M.A. Luersen, R. Le Riche, Globalized Nelder-Mead method for engineering optimization, Comput. Struct. 82 (2004) 2251–2260. [19] R. Barati, Parameter estimation of nonlinear Muskingum models using Nelder-Mead simplex algorithm, J. Hydrol. Eng. 16 (2011) 946–954. [20] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313. [21] A. Ouria, M.M. Toufigh, Application of Nelder-Mead simplex method for unconfined seepage problems, Appl. Math. Model. 33 (2009) 3589–3598. [22] X.J. Tan, W.Z. Chen, H.M. Tian, et al., Water flow and heat transport including ice/ water phase change in porous media: numerical simulation and application, Cold Reg. Sci. Technol. 68 (2011) 74–84. [23] D.R. Lide, Handbook of Chemistry and Physics, 84th ed., CRC Press Inc., 2003–2004, pp. 980–981. [24] X.J. Tan, W.Z. Chen, D.S. Yang, et al., Study on the influence of airflow on the temperature of the surrounding rock in a cold region tunnel and its application to insulation layer design, Appl. Therm. Eng. 67 (2014) 320–334.

It can be seen that our combining simulation method is suitable for circular section excavation, around which the freeze pipes are also arranged in a circle. However, there are many shapes of cavern including U-shaped, arch, square section and so on [3,34]. The optimum distribution of freeze pipes are more complicate in those cases. Therefore, it’s necessary to establish the possible distribution function of optimal positions of freeze pipes around non-circular cross section. (2) Combining optimization of artificial ground freezing with other presupported technology It has been derived the freezing time exceeds 38 days after optimization of freeze pipes arrangement when the initial seepage flow is more than 2.0 m/d in the case study as is shown in Fig. 14. Therefore, the proposed optimization method combining with other pre-supported method may be more economical and time-saving, such as the combining method of pipe-roofing and artificial ground freezing [14]. The application of tubular roof between the freeze pipes could well prevent water flow, enhance heat transfer and decrease the deformation of surrounding environments [15]. Therefore, a more effective and energy-saving optimization method of artificial ground freezing combined with other pre-supported method may need to be studied further. 6. Summary A new combining optimization method of freeze pipes arrangement considering water flow has been proposed and applied to the artificial ground freezing process during tunnel excavation. From this study, the following conclusions can be drawn: (1) The developed coupled hydro-thermal model considering water/ice phase transition is suitable to simulate the freezing process of porous media during artificial ground freezing under water flow condition. (2) Through the comparisons between the calculated and measured temperature, the unfrozen water content and the effective thermal conductivity of freezing porous media could be well expressed by the exponential weighted mean model and the exponential function, respectively. (3) The development of freezing band around freeze pipes is significantly influenced by water flow and the completed freezing time of designed arch tremendously increases with the increase in the initial flow velocity of water, which can be well predicted by an exponential equation. (4) The proposed combining optimization method is very appropriate for searching out the most optimum positions of freeze pipes under seepage flow no matter which one of the common distribution 444

Applied Thermal Engineering 135 (2018) 435–445

S. Huang et al.

[32] E. Pimentel, A. Sres, G. Anagnostou, Large-scale laboratory tests on artificial ground freezing under seepage-flow conditions, Geotechnique 62 (2012) 227–241. [33] J.M. Konrad, Frost susceptibility related to soil index properties, Can. Geotech. J. 36 (1999) 403–417. [34] S.Y. Li, Y.M. Lai, M.Y. Zhang, et al., Minimum ground pre-freezing time before excavation of Guangzhou subway tunnel, Cold Reg. Sci. Technol. 46 (2006) 181–191. [35] Y. Yamamoto, S.M. Springman, Axial compression stress path tests on artificial frozen soil samples in a triaxial device at temperatures just below 0°C, Can. Geotech. J. 51 (2014) 1178–1195. [36] X.D. Hu, J.Z. Yu, H. Ren, et al., Analytical solution to steady-state temperature field for straight-row-piped freezing based on superposition of thermal potential, Appl. Therm. Eng. 111 (2017) 223–231. [37] R. Lackner, A. Amon, H. Lagger, Artificial ground freezing of fully saturated soil: thermal problem, J. Eng. Mech. 131 (2005) 211–220.

[25] S.B. Huang, Q.S. Liu, A.P. Cheng, et al., A fully coupled thermo-hydro-mechanical model including the determination of coupling parameters for freezing rock, Int. J. Rock Mech. Min. 103 (2018) 205–2144. [26] S.B. Huang, Q.S. Liu, Y.Z. Liu, et al., Freezing strain model for estimating the unfrozen water content of saturated rock under low temperature, Int. J. Geomech. 18 (2018) 04017137, http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001057. [27] J.M. Konrad, N.R. Morgenstern, The segregation potential of a freezing soil, Can. Geotech. J. 18 (1981) 482–491. [28] J.M. Konrad, Influence of cooling rate on the temperature of ice lens formation in clayey silts, Cold Reg. Sci. Technol. 16 (1989) 25–36. [29] J.M. Konrad, N.R. Morgenstern, Effects of applied pressure on freezing soils, Can. Geotech. J. 19 (1982) 494–505. [30] J.M. Konrad, M. Shen, 2-D frost action modeling using the segregation potential of soils, Cold Reg. Sci. Technol. 24 (1996) 263–278. [31] S. Nishimura, A. Gens, S. Olivella, et al., THM-coupled finite element analysis of frozen soil: formulation and application, Geotechnique 59 (2009) 159–171.

445