# Thermo-Hydromechanical Coupling Responses Driven by a Central Heat Source in Unsaturated Soils

Authors: Wei Wang; Jingjing Liu; Jing Chen
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

Based on the thermodynamic and thermoelastic theory, the coupled governing equations of deformation, heat transfer, and moisture migration in unsaturated soils were given. The coupled calculation process was realized by the finite element method. The hydrothermal coupling characteristics of two cases of unsaturated kaolin clay were studied by using a self-developed test device, and the test results were compared with the numerical results. The results showed that when the initial saturation of the soil is high, the volumetric water content of the measured point increases, and then it is in a stable state during the heating process. When the initial saturation is low, the volumetric water content near the heat source increases firstly and then decreases during the heating process. The rise and fall of temperature make the volumetric water content of the soil irreversible. The volumetric water content of each measured point is lower than the initial state. The closer it is to the internal measured heat source, the more obvious this phenomenon is; the lower the initial saturation is, the more obvious it is.


---

## Body

## 1. Introduction

The study on the unsteady responses of internal temperature transfer and moisture migration induced by the coupling effect of suction potential, gravitational potential, and temperature potential in unsaturated soils has very important significance. The application respects are involved in geothermal resource development [1], thermal energy storage [2], nuclear waste disposal [3], and heating pipeline design [4]. The study mainly relates to the coupling of different physical fields, for instance, hydraulic and mechanical thermal responses.In order to analyse the responses of various physical fields, scholars have studied them from different angles. Chen et al. [5] constructed a complicated multifield coupling model for unsaturated porous media considering chemical and hydration swelling effects based on mixture coupling theory, nonequilibrium dynamics, and Biot’s elastic theory. Chen et al. [6] proposed a method for analysing thermo-hydromechanical processes in porous media based on the continuum mechanics of the conservation of momentum, energy, and mass and the averaging method in the three-phase mixture theory (solid-liquid-gas) system. This model considered the effects of stress-strain, water-air-vapour phase flow, thermal transfer, and porosity evolution on the coupled response. Rutqvist et al. [7] gave basic governing equations of the thermo-hydromechanical coupling process of unsaturated fractured rock and soil and compared four finite element programs. Chen et al. [8] conducted experimental evaluations on several models describing the relationship between permeability and capillary pressure-saturation, which showed that LNM, VGM, GDM, and BRB models can better describe the flow of two-phase fluids. Wu et al. [9] established a coupled thermo-hydromechanical stress-strain relationship of unsaturated soils, which can describe the influence of heat on the hydraulic characteristics of soils and can quantitatively simulate thermal softening phenomena. Finally, the test results were compared with the numerical simulation results. Akrouch et al. [10] discussed the coupled hydrothermal characteristics of energy piles to unsaturated soils by experiments, theory, and numerical simulation and compared the experimental results with analytical and finite element solutions. Ghorbani et al. [11] established a theory that can describe the hydraulic interaction behavior for typical unsaturated soils under static and dynamic loads on the basis of previous studies. By solving the dynamic problems of unsaturated soils with elastic frameworks and nonlinear frameworks, the static problem of the saturated soil verified the correctness of the established model.Bai [12] discussed the semianalytical solutions of porous and saturated spheres with cavities under several boundary conditions and compared the effects of thermal response on field variables such as radial displacement, pore pressure, and tangential or radial stress through numerical analysis, which is a pioneering work. Caulk et al. [13] proposed a pore-scale numerical method for simulating the thermo-hydromechanical coupling behavior of granular media on the basis of the existing fluid mechanics model. For unsaturated soils, Takayama et al. [14] established a constitutive model based on the Cambridge model to describe the changes of the mechanical properties in view of bentonite in the nuclear waste repository due to the hydration process. Zeng et al. [15] established a vertical one-dimensional, mass flow heat, and two-phase model considering diffusion, advection, and dispersion and used the model to discuss the advection effects of high-permeability soil as well as low-permeability soil. For saturated porous media, Bai [16] established a normalized coupled governing equation and discussed the theoretical response under cyclic thermal load for one-dimensional cyclic thermal consolidation. Among these, the evolution process of displacement, pore pressure, and temperature from an unsteady state (i.e., instantaneous state) to a quasi-steady state was analysed. An et al. [17] established a two-dimensional hydrothermal coupling model and used meteorological data to estimate the surface heat flow and water flow boundary. On this basis, the variation rules of temperature as well as the volumetric water content of embankment were also studied. Zhou and Ng [18] established a new coupled thermomechanical theory of saturated soil based on the boundary surface plasticity theory. The shape and size of the boundary surface of this theoretical model can change with temperature, which can predict the elastic-plastic response of saturated soil under small deformation. For unsaturated porous media, the multifield coupling process is greatly affected by saturation degree. The effect of saturation on the coupled thermo-hydromechanical process in porous media is of great significance to understand the properties of the multifield coupling response, even including the influence of chemical factors [19–21]. However, all the constitutive models described above are derived based on the general elastic-plastic theory.In this paper, the relationship between the degree of saturation and pore gas pressure as well as pore water pressure in the soil-water characteristic curve [22] is utilized to derive the seepage equation of unsaturated porous media through the continuity equation of pore water and pore gas. Based on the combination of thermodynamics and thermoelastic theory, the governing equations with coupled thermo-hydromechanical characteristics as well as state equations of the water-gas two-phase flow in porous media are further given. The coupling calculation process is realized by the finite element method. Then, a self-developed test device was utilized to analyse the hydrothermal coupling characteristics of two different unsaturated kaolin clay under the same temperature load, and the obtained experimental results verified the rationality of the proposed theoretical model. This work provides further insight into the hydrothermal effects of kaolin soil from the multifield coupling point of view.


## 2. Coupled Governing Equations for Unsaturated Soils

### 2.1. Seepage Differential Equation

Under the condition of ignoring the evaporation of water in porous media, the mass conservation equation of pore water, pore gas, and solid matrix can be expressed as [7, 8](1)∂nSwρw∂t+∇⋅nSwρw⋅vw=0,(2)∂nSgρg∂t+∇⋅nSgρg⋅vg=0,(3)∂1−nρs∂t+∇⋅1−nρs⋅vs=0,where ∇·is the divergence calculation; n is the porosity of the porous medium; ρw, ρs, and ρg are the density of pore water, solid matrix, and pore gas (ML−3), respectively; Sw and Sg are pore water and pore gas saturation; and vw, vs, and vg are the velocity vectors of pore water, solid substrate, and pore gas (LT−1), respectively.Considering the effect of gravity, the velocity of the pore water as well as pore gas phase to the solid substrate can be expressed as [23–26](4)vrw=vw−vs,vrg=vg−vs,∇⋅vs≈∂εv∂t.According to [27–31],(5)nSwvrw=−K⋅Krwμw⋅∇Pw+ρwg⋅∇z=uw,(6)nSgvrg=−K⋅Krgμg⋅∇Pg+ρgg⋅∇z=ug,where K is the porous medium permeability (L2); μw and μg are the pore water and pore gas viscosity (ML−1T−1), respectively; Krw and Krg are the relative permeability of the liquid and gas phase, respectively; Pw and Pg are the pore water and pore gas pressure (ML−1T−2); vrw and vrg are the relative velocity of the liquid and gas phase, respectively; z is the height relative to the reference plane; g is the gravitational acceleration; ug and uw are the average relative flow rates of pore gas and pore water (LT−1) in the cross section of the porous medium, respectively.According to [32–36], it can be known that pore water saturation should be a function of temperature and matrix suction, which has the following relationship:(7)Sw=SwPc,T,Sg=1−Sw,Pc=Pg−Pw,where T is the current absolute temperature and Pc is the matrix suction (ML−1T−2).Then, the derivative of the saturation of pore waterSw with respect to time is(8)∂Sw∂t=∂Sw∂Pc∂Pc∂t+∂Sw∂T∂T∂t=CP∂Pg∂t−CP∂Pw∂t+CT∂T∂t=−∂Sg∂t,where CP and CT are the water-holding characteristic coefficients of capillary action (M−1LT2) and temperature action (θ−1), respectively, which are determined by the soil-water characteristic curve.Substituting equations (4)–(8) into equations (1) and (2), the seepage differential equations of pore water as well as pore gas can be obtained as follows:(9)nρwCP∂Pw∂t−∇⋅−ρwKKrwμw∇Pw+ρwg⋅∇z=nρwCP∂Pg∂t+nρwCT∂T∂t+nSw∂ρw∂t+Swρw1−nρs∂ρs∂t+∂εv∂t,nρgCP∂Pg∂t−∇⋅−ρgKKrgμg∇Pg+ρgg⋅∇z=nρgCP∂Pw∂t−nρgCT∂T∂t+nSg∂ρg∂t+Sgρg1−nρs∂ρs∂t+∂εv∂t.


### 2.2. Momentum Conservation Equation

For unsaturated porous media, the momentum conservation equation without considering the inertial force can be expressed as [7, 8](10)∇⋅σ+1−nρs+nSwρw+nSgρg⋅g=0.Considering the skeleton deformation, thermal strain, wet strain, pore pressure deformation, and equivalent stress principles of porous media, the stress-strain relationship is(11)σ=D:ε−εT−εS−εP−IP¯,where(12)ε=TβTT−T0,ε=SβMSw−Sw0,εP=−P¯−P¯03Km,P¯=SwPw+SgPg,where βT and βM are the linear thermal expansion coefficient and linear wet expansion coefficient of the porous solid matrix, respectively; D is the tangential stiffness matrix, which is the weighted average pore pressure (ML−1T−2); Sw0 and T0 are the initial saturation and initial temperature, respectively; and Km is the bulk modulus of the matrix of the porous medium (ML−1T−2).


### 2.3. Energy Conservation Equation

Considering the change rate of thermal strain and wet strain energy, for unsaturated porous media, the energy conservation equation can be expressed as [7, 8](13)∂∂tnSwρwew+nSgρgeg+1−nρses+1−n3KmβTT+βMSw∂εv∂t+∇⋅kt⋅∇T+nSwPw∇⋅vw+nSgPg∇⋅vg+∇⋅hwnSwρw⋅vw+hgnSgρg⋅vg+hs1−nρs⋅vs=Q,where ew, eg, and es are the internal energy per unit mass of pore water, pore gas, and solid matrix (L2T−2), respectively, e=cv·(T − T0); hw, hg, and hs are the heat content of pore water, pore gas, and solid matrix per unit mass, respectively (L2T−2); Q is the total heat source intensity (L−1MT−3); kt is the total thermal conductivity of porous media (MT−3θ−1); cv and cp are the specific heat of constant volume and pressure, respectively (L2T−2θ−1); and γ is the specific heat rate.On the left side of equation (13), the first term is the change rate of the internal energy of the system, the second term is the change rate of reversible elastic energy caused by the solid phase change [37–40], the third term is the heat conduction, the fourth term is the energy change produced by pore water pressure [41–44], and the fifth term is the difference between the incoming and outgoing energy per unit time.According to previous studies [45–49], considering the effects of saturation and porosity degree on the effectivity of the thermal conductivity, the total thermal conductivity can be expressed as(14)kt=na1kw−kgSw1+a2Sw+kg+1−nks,where ks, kg, and kw are the thermal conductivity of the porous media matrix, pore gas, and pore water (MT−3θ−1) and a1 and a2 both are experimental fitting parameters, respectively.


### 2.4. Water Retention Characteristics

In view of the VG model [50–53], the relationship between matrix suction and relative permeability and relative effective saturation can be expressed as(15)Sew=11+αHCnmHC>0,1HC≤0,HC=Pcρwg,Krw=SewL1−1−Sew1/mm2HC>0,0HC≤0,Krg=1−SewL1−Sew1/m2mHC>0,1HC≤0,where Sew is the relative effective saturation; m and n are the fitting parameters, where m is related to the asymmetry of the soil-water characteristic curve and n is related to the slope of the curve; and Hc is the matrix suction head.Pore water saturation degree can be written as(16)Sw=Swr+Sws−SwrSew,where Swr and Sws are the residual saturation and maximum saturation, respectively.According to equation (8), it is assumed that the parameters α, m, and L have approximately linear relationships with temperature, which can be commonly written as(17)αα0=1+CαTT−T0,mm0=1+CmTT−T0,LL0=1+CLTT−T0.WhenHc > 0, according to equations (8), (14), (15), and (16), one has(18)CP=∂Sw∂t=−αmSws−Swrρwg1−mSew1/m1−Sew1/mm,CT=∂Sw∂T=CmTm0SwSws−Swr1−m1m−Sw1/mlnSw−m1−Sw1/mln1−Sw1/m−CαTα0mSwSws−Swrα1−m1−Sw1/m.


### 2.5. State Equations

In order to simulate the real situation effectively, each state quantity adopts a nonlinear calculation formula related to temperature. According to [54–56], the density ρ of pore gas, pore water, and solid matrix, the viscosity coefficient of pore gas and pore water μ (Pa·s), the heat capacity at constant pressure cp (J/(kg·K)), and thermal conductivity k (W/(m·K)) can be expressed as, respectively,(19)ρw=838.47+1.40T−0.003T3+3.72×10−7T3,ρg=PgRgT,ρs=ρs01+P¯−P¯0Km−3βTT−T0−3βMSw−Sw0−trσ′−σ0′3Km1−n,μw=1.38−0.02T+1.36T2−4.65×10−7T3+8.90×10−10T4−9.08×10−13T5+3.85×10−16T6273.15K<T<413.15K,μg=−8.38×10−7+8.36×10−8T−7.69×10−11T2+4.64×10−14T3−1.07×10−17T4203.15K<T<1573.15K,cpw=12010.15−80.41T+0.31T2−5.38×10−4T3+3.63×10−7T4273.15K<T<553.15K,cpg=1047.64−0.37T+9.45×10−4T2−6.02×10−7T3+1.29×10−10T4203.15K<T<1573.15K,kw=−0.87+8.95×10−3T−1.58×10−5T2−7.98×10−9T3273.15K<T<553.15K,kg=−2.28×10−3+1.16×10−4T−7.90×10−8T2−4.12×10−11T3−7.44×10−15T4203.15K<T<1573.15K,where(20)trσ′−σ0′=3Kbεv+3KbKmP¯−P¯0−9KbβTT−T0−9KbβMSw−Sw0,where Rg is actually the gas constant of air (L2T−2θ−1); σ′ is the effective stress acting on the matrix in the porous medium (L−1MT−2), in which σ0′ is the initial value of the effective stress (L−1MT−2); ρs0 is the initial density of the porous medium matrix (ML−3); and tr(⋅) is the trace operation in mathematics.According to equations (3), (11), and (19), the porosity can be expressed as(21)n=1−1−n0eεv0−εv1+P¯−P¯0/Km−3βTT−T0−3βMSw−Sw0−trσ′−σ0′/3Km1−n.


## 2.1. Seepage Differential Equation

Under the condition of ignoring the evaporation of water in porous media, the mass conservation equation of pore water, pore gas, and solid matrix can be expressed as [7, 8](1)∂nSwρw∂t+∇⋅nSwρw⋅vw=0,(2)∂nSgρg∂t+∇⋅nSgρg⋅vg=0,(3)∂1−nρs∂t+∇⋅1−nρs⋅vs=0,where ∇·is the divergence calculation; n is the porosity of the porous medium; ρw, ρs, and ρg are the density of pore water, solid matrix, and pore gas (ML−3), respectively; Sw and Sg are pore water and pore gas saturation; and vw, vs, and vg are the velocity vectors of pore water, solid substrate, and pore gas (LT−1), respectively.Considering the effect of gravity, the velocity of the pore water as well as pore gas phase to the solid substrate can be expressed as [23–26](4)vrw=vw−vs,vrg=vg−vs,∇⋅vs≈∂εv∂t.According to [27–31],(5)nSwvrw=−K⋅Krwμw⋅∇Pw+ρwg⋅∇z=uw,(6)nSgvrg=−K⋅Krgμg⋅∇Pg+ρgg⋅∇z=ug,where K is the porous medium permeability (L2); μw and μg are the pore water and pore gas viscosity (ML−1T−1), respectively; Krw and Krg are the relative permeability of the liquid and gas phase, respectively; Pw and Pg are the pore water and pore gas pressure (ML−1T−2); vrw and vrg are the relative velocity of the liquid and gas phase, respectively; z is the height relative to the reference plane; g is the gravitational acceleration; ug and uw are the average relative flow rates of pore gas and pore water (LT−1) in the cross section of the porous medium, respectively.According to [32–36], it can be known that pore water saturation should be a function of temperature and matrix suction, which has the following relationship:(7)Sw=SwPc,T,Sg=1−Sw,Pc=Pg−Pw,where T is the current absolute temperature and Pc is the matrix suction (ML−1T−2).Then, the derivative of the saturation of pore waterSw with respect to time is(8)∂Sw∂t=∂Sw∂Pc∂Pc∂t+∂Sw∂T∂T∂t=CP∂Pg∂t−CP∂Pw∂t+CT∂T∂t=−∂Sg∂t,where CP and CT are the water-holding characteristic coefficients of capillary action (M−1LT2) and temperature action (θ−1), respectively, which are determined by the soil-water characteristic curve.Substituting equations (4)–(8) into equations (1) and (2), the seepage differential equations of pore water as well as pore gas can be obtained as follows:(9)nρwCP∂Pw∂t−∇⋅−ρwKKrwμw∇Pw+ρwg⋅∇z=nρwCP∂Pg∂t+nρwCT∂T∂t+nSw∂ρw∂t+Swρw1−nρs∂ρs∂t+∂εv∂t,nρgCP∂Pg∂t−∇⋅−ρgKKrgμg∇Pg+ρgg⋅∇z=nρgCP∂Pw∂t−nρgCT∂T∂t+nSg∂ρg∂t+Sgρg1−nρs∂ρs∂t+∂εv∂t.


## 2.2. Momentum Conservation Equation

For unsaturated porous media, the momentum conservation equation without considering the inertial force can be expressed as [7, 8](10)∇⋅σ+1−nρs+nSwρw+nSgρg⋅g=0.Considering the skeleton deformation, thermal strain, wet strain, pore pressure deformation, and equivalent stress principles of porous media, the stress-strain relationship is(11)σ=D:ε−εT−εS−εP−IP¯,where(12)ε=TβTT−T0,ε=SβMSw−Sw0,εP=−P¯−P¯03Km,P¯=SwPw+SgPg,where βT and βM are the linear thermal expansion coefficient and linear wet expansion coefficient of the porous solid matrix, respectively; D is the tangential stiffness matrix, which is the weighted average pore pressure (ML−1T−2); Sw0 and T0 are the initial saturation and initial temperature, respectively; and Km is the bulk modulus of the matrix of the porous medium (ML−1T−2).


## 2.3. Energy Conservation Equation

Considering the change rate of thermal strain and wet strain energy, for unsaturated porous media, the energy conservation equation can be expressed as [7, 8](13)∂∂tnSwρwew+nSgρgeg+1−nρses+1−n3KmβTT+βMSw∂εv∂t+∇⋅kt⋅∇T+nSwPw∇⋅vw+nSgPg∇⋅vg+∇⋅hwnSwρw⋅vw+hgnSgρg⋅vg+hs1−nρs⋅vs=Q,where ew, eg, and es are the internal energy per unit mass of pore water, pore gas, and solid matrix (L2T−2), respectively, e=cv·(T − T0); hw, hg, and hs are the heat content of pore water, pore gas, and solid matrix per unit mass, respectively (L2T−2); Q is the total heat source intensity (L−1MT−3); kt is the total thermal conductivity of porous media (MT−3θ−1); cv and cp are the specific heat of constant volume and pressure, respectively (L2T−2θ−1); and γ is the specific heat rate.On the left side of equation (13), the first term is the change rate of the internal energy of the system, the second term is the change rate of reversible elastic energy caused by the solid phase change [37–40], the third term is the heat conduction, the fourth term is the energy change produced by pore water pressure [41–44], and the fifth term is the difference between the incoming and outgoing energy per unit time.According to previous studies [45–49], considering the effects of saturation and porosity degree on the effectivity of the thermal conductivity, the total thermal conductivity can be expressed as(14)kt=na1kw−kgSw1+a2Sw+kg+1−nks,where ks, kg, and kw are the thermal conductivity of the porous media matrix, pore gas, and pore water (MT−3θ−1) and a1 and a2 both are experimental fitting parameters, respectively.


## 2.4. Water Retention Characteristics

In view of the VG model [50–53], the relationship between matrix suction and relative permeability and relative effective saturation can be expressed as(15)Sew=11+αHCnmHC>0,1HC≤0,HC=Pcρwg,Krw=SewL1−1−Sew1/mm2HC>0,0HC≤0,Krg=1−SewL1−Sew1/m2mHC>0,1HC≤0,where Sew is the relative effective saturation; m and n are the fitting parameters, where m is related to the asymmetry of the soil-water characteristic curve and n is related to the slope of the curve; and Hc is the matrix suction head.Pore water saturation degree can be written as(16)Sw=Swr+Sws−SwrSew,where Swr and Sws are the residual saturation and maximum saturation, respectively.According to equation (8), it is assumed that the parameters α, m, and L have approximately linear relationships with temperature, which can be commonly written as(17)αα0=1+CαTT−T0,mm0=1+CmTT−T0,LL0=1+CLTT−T0.WhenHc > 0, according to equations (8), (14), (15), and (16), one has(18)CP=∂Sw∂t=−αmSws−Swrρwg1−mSew1/m1−Sew1/mm,CT=∂Sw∂T=CmTm0SwSws−Swr1−m1m−Sw1/mlnSw−m1−Sw1/mln1−Sw1/m−CαTα0mSwSws−Swrα1−m1−Sw1/m.


## 2.5. State Equations

In order to simulate the real situation effectively, each state quantity adopts a nonlinear calculation formula related to temperature. According to [54–56], the density ρ of pore gas, pore water, and solid matrix, the viscosity coefficient of pore gas and pore water μ (Pa·s), the heat capacity at constant pressure cp (J/(kg·K)), and thermal conductivity k (W/(m·K)) can be expressed as, respectively,(19)ρw=838.47+1.40T−0.003T3+3.72×10−7T3,ρg=PgRgT,ρs=ρs01+P¯−P¯0Km−3βTT−T0−3βMSw−Sw0−trσ′−σ0′3Km1−n,μw=1.38−0.02T+1.36T2−4.65×10−7T3+8.90×10−10T4−9.08×10−13T5+3.85×10−16T6273.15K<T<413.15K,μg=−8.38×10−7+8.36×10−8T−7.69×10−11T2+4.64×10−14T3−1.07×10−17T4203.15K<T<1573.15K,cpw=12010.15−80.41T+0.31T2−5.38×10−4T3+3.63×10−7T4273.15K<T<553.15K,cpg=1047.64−0.37T+9.45×10−4T2−6.02×10−7T3+1.29×10−10T4203.15K<T<1573.15K,kw=−0.87+8.95×10−3T−1.58×10−5T2−7.98×10−9T3273.15K<T<553.15K,kg=−2.28×10−3+1.16×10−4T−7.90×10−8T2−4.12×10−11T3−7.44×10−15T4203.15K<T<1573.15K,where(20)trσ′−σ0′=3Kbεv+3KbKmP¯−P¯0−9KbβTT−T0−9KbβMSw−Sw0,where Rg is actually the gas constant of air (L2T−2θ−1); σ′ is the effective stress acting on the matrix in the porous medium (L−1MT−2), in which σ0′ is the initial value of the effective stress (L−1MT−2); ρs0 is the initial density of the porous medium matrix (ML−3); and tr(⋅) is the trace operation in mathematics.According to equations (3), (11), and (19), the porosity can be expressed as(21)n=1−1−n0eεv0−εv1+P¯−P¯0/Km−3βTT−T0−3βMSw−Sw0−trσ′−σ0′/3Km1−n.


## 3. Laboratory Model Experiment

### 3.1. Test Apparatus Description

A test device with a cylindrical chamber is used to study the water-heat coupling effect in unsaturated soils, as shown in Figure1. The upper and lower ends of the test device are fixedly connected with the stainless steel plate through movable bolts, and an O-shaped rubber ring is set for sealing treatment. The test device is made of stainless steel (304 stainless steel), with an inner radius of R = 50 cm, a cylinder of 1 cm wall thickness, and a height of H = 50 cm. In order to meet the axisymmetric temperature, hydraulic, deformation, and migration boundary conditions, the stainless steel bottom plate and top cover are inlaid with a heat-insulating and impervious PC board. A layer of asbestos material is placed between the top cover and the bottom plate and the PC board. The centre of the test device is a cylindrical heat source Tin with a radius of r = 2 cm, whose temperature is applied by using a constant temperature water bath circulation device, which can provide a long-term constant temperature possessing an accuracy of ±0.1°C. Besides, the outer temperature boundary of the test device also circulates through the constant temperature water bath. The device is applied and set to Tout = 25°C.Figure 1
Schematic diagram of the test device.In order to analyse the heat-water coupling effect in unsaturated soils, temperature sensors (K-type thermocouples) and volumetric water content sensors (EC-5) were embedded at different radial distances on the middle plane of the soil during the filling process of the soil sample, as shown in Figure 1. Four measuring points of temperature and moisture were arranged from the edge of the test tube to the centre, with an interval of 10 cm. By the temperature and moisture collector, the moisture and temperature changes at different locations were recorded in real time.


### 3.2. Test Schemes

In the experiment, washed kaolin clay from Lingshou County, Hebei Province, China, was used as the porous medium, and the relevant property parameters are shown in Table1. A mortar mixer is used to mix kaolin clay with a particle size of 37.5 μm into two different saturation degrees. Through experiments, the changes in temperature and seepage field of kaolin clay with two different saturation degrees are studied during the heating and cooling process. The two cases include (1) the kaolin sample with the initial saturation of 0.47 (mass moisture content 0.130), the internal temperature boundary of 80°C, and the outer temperature boundary of 25°C. The coupled properties of the seepage and temperature field during the heating process are studied. Both internal and external temperature boundaries are at 25°C. Meanwhile, the coupling response of seepage and temperature field is considered during the cooling process; (2) the initial saturation of the kaolin sample is 0.29 (mass moisture content 0.08), the internally measured temperature boundary is 80°C, and the outer temperature boundary is 25°C. The schemes are emerged in Table 2.Table 1
Related physical parameters of kaolin.
Processing typeDensity (g/cm3)ColourRefractorinessSiO2 content (%)Molecular weightWater wash2.56White170052258.16Table 2
Test schemes.
Serial numberInitial saturationInitial volumetric water contentTin (°C)Tout (°C)Dry density (g/cm3)10.470.19680 (heating)25 (cooling)251.520.290.12080 (heating)25 (cooling)251.5


## 3.1. Test Apparatus Description

A test device with a cylindrical chamber is used to study the water-heat coupling effect in unsaturated soils, as shown in Figure1. The upper and lower ends of the test device are fixedly connected with the stainless steel plate through movable bolts, and an O-shaped rubber ring is set for sealing treatment. The test device is made of stainless steel (304 stainless steel), with an inner radius of R = 50 cm, a cylinder of 1 cm wall thickness, and a height of H = 50 cm. In order to meet the axisymmetric temperature, hydraulic, deformation, and migration boundary conditions, the stainless steel bottom plate and top cover are inlaid with a heat-insulating and impervious PC board. A layer of asbestos material is placed between the top cover and the bottom plate and the PC board. The centre of the test device is a cylindrical heat source Tin with a radius of r = 2 cm, whose temperature is applied by using a constant temperature water bath circulation device, which can provide a long-term constant temperature possessing an accuracy of ±0.1°C. Besides, the outer temperature boundary of the test device also circulates through the constant temperature water bath. The device is applied and set to Tout = 25°C.Figure 1
Schematic diagram of the test device.In order to analyse the heat-water coupling effect in unsaturated soils, temperature sensors (K-type thermocouples) and volumetric water content sensors (EC-5) were embedded at different radial distances on the middle plane of the soil during the filling process of the soil sample, as shown in Figure 1. Four measuring points of temperature and moisture were arranged from the edge of the test tube to the centre, with an interval of 10 cm. By the temperature and moisture collector, the moisture and temperature changes at different locations were recorded in real time.


## 3.2. Test Schemes

In the experiment, washed kaolin clay from Lingshou County, Hebei Province, China, was used as the porous medium, and the relevant property parameters are shown in Table1. A mortar mixer is used to mix kaolin clay with a particle size of 37.5 μm into two different saturation degrees. Through experiments, the changes in temperature and seepage field of kaolin clay with two different saturation degrees are studied during the heating and cooling process. The two cases include (1) the kaolin sample with the initial saturation of 0.47 (mass moisture content 0.130), the internal temperature boundary of 80°C, and the outer temperature boundary of 25°C. The coupled properties of the seepage and temperature field during the heating process are studied. Both internal and external temperature boundaries are at 25°C. Meanwhile, the coupling response of seepage and temperature field is considered during the cooling process; (2) the initial saturation of the kaolin sample is 0.29 (mass moisture content 0.08), the internally measured temperature boundary is 80°C, and the outer temperature boundary is 25°C. The schemes are emerged in Table 2.Table 1
Related physical parameters of kaolin.
Processing typeDensity (g/cm3)ColourRefractorinessSiO2 content (%)Molecular weightWater wash2.56White170052258.16Table 2
Test schemes.
Serial numberInitial saturationInitial volumetric water contentTin (°C)Tout (°C)Dry density (g/cm3)10.470.19680 (heating)25 (cooling)251.520.290.12080 (heating)25 (cooling)251.5


## 4. Results and Discussion

### 4.1. Physical Model and Boundary Conditions

This paper uses “COMSOL Multiphysics 5.6” finite element software to solve the governing equations in Section2 and to simulate the experimental process. The calculation model adopts an axisymmetric model. The inner heat source and the outer cylindrical wall were set as impervious and incapable of displacement boundary conditions; the top cover and bottom plate were set as thermally insulated with impervious boundary conditions. The temperature boundary of the inner heat source was set to 80°C and 25°C during the rising and falling process. The temperature boundary of the outer cylindrical wall was set to 25°C, and the initial temperature of the soil was 25°C. Thus, the optimization results of relevant parameters required in the calculation are shown in Table 3. The model mesh was divided into a free triangle mesh, the maximum element size was 2.65 cm, and the minimum element size was 0.15 mm (Figure 2). The complete grid contains a total of 884 domain elements and 74 boundary elements.Table 3
Model parameters.
Parameterρs0Ea1a2ksKcpsβTβMg/cm3MPaW/(m·K)m2J/(kg·K)1/KValue1.518.321.220.681.0 × 10−229702 × 10−9−1×10−3Figure 2
The meshing of the computational model.In the calculation process, both the internal heat source temperature boundaryTin and the outer temperature boundary Tout were set as the first type of boundary conditions. As a result, the internal heat source temperature boundary is within the first 5.15 d, and the temperature was maintained at about 80°C. During the test, the temperature was basically maintained at about 25°C. It should be noted that, during the 5.15 d cooling process, the water bath needs to be replaced with cold water to achieve the cooling effect. The constant temperature water bath had no cooling function, so the internal temperature of the heat source was measured. After the temperature drop starts, it first rises and then gradually decreases to 25°C, as shown in Figure 3(a), while the outer temperature boundary has been maintained at about 25°C (Figure 3(b)), which was arranged on the outer cylindrical wall of the test device (Figure 3(b)). The temperature measured by sensors changes over time at four different locations.Figure 3
Variation of (a) heat source temperature and (b) boundary temperature with time.
(a)(b)


### 4.2. Discussion

According to the test results, it can be found that the volumetric water content of each monitoring position in the soil has changed after the temperature rising and falling process. Due to different soil-water characteristics of unsaturated soils in the process of moisture absorption and dedrying, the VG model parameters are different during the heating and cooling processes. After trial calculations, the parameterα1 = 1.89 m−1 during the heating process; it is taken as α2 = 1.93 m−1 during the cooling process. In the calculation process, the initial pore water pressure as well as pore gas pressure was 0.94 × 105 Pa and 1.01 × 105 Pa, respectively. Figure 3 shows the comparison between the calculated results and the experimental results of the volumetric water content at different monitoring positions in the soil for case 1. By comparison, it can be observed that the calculation results are consistent with the experimental results. During the heating process, the calculation results of positions B1, C1, and D1 are in good agreement with the test results, while the calculation results of position A1 are slightly lower than the test results, but the trend of change is still consistent. It can be observed that when the parameter α does not change during the heating and cooling process (α1 = α2 = 1.89 m−1) by comparing Figures 4(a) and 4(b), the calculation results present that the content of volumetric water at each location after the heating and cooling process has little change from the initial state and is far from the test result (Figure 4(a)). When the parameter α is set differently during the heating and cooling process (e.g., α1 = 1.89 m−1 and α2 = 1.93 m−1), the calculation results show that the volumetric water content at each location after the heating and cooling process is lower than the initial state and closer to the test results (Figure 4(b)).Figure 4
Comparison of the volumetric water content with numerical calculation in case 1: (a)α1 = α2 = 1.89 m−1; (b) α1 = 1.89 m−1 and α2 = 1.93 m−1.
(a)(b)Figure5 shows the comparison between the test results and the calculated results of the temperature at different monitoring positions in the soil for cases 1 and 2. In the calculation process, the parameter α of case 1 is set to α1 = 1.89 m−1 and α2 = 1.93 m−1, and the parameter α of case 2 is set to α1 = 1.89 m−1 and α2 = 1.94 m−1. By comparison, it can be seen that the numerical calculation results of the temperature changes at each measuring point during the heating process of the two cases are in good agreement with the test results, but there is a certain gap during the cooling process. This is due to the errors in the parameters and boundary conditions of the laboratory experiments and numerical calculations.Figure 5
Comparison of results of temperature with numerical calculations in case 1: (a)α1 = 1.89 m−1 and α2 = 1.93 m−1; (b) α1 = 1.89 m−1 and α2 = 1.94 m−1.
(a)(b)Figure6 shows the comparison between the test results and the calculated results of the volumetric water content at different monitoring positions in the soil for case 2. Through comparison, it can be observed that the calculation results are consistent with the experimental results. In the calculation process, the parameter α is set to α1 = 1.89 m−1 and α2 = 1.94 m−1, and the initial pore water pressure and pore gas pressure are set to 1.56 × 105 Pa and 1.46 × 105 Pa, respectively. The permeability K in saturated soil is set to 4 × 10−19 m2. It can be observed from Figure 6 that the numerical calculation results in the heating process are in good agreement with the test results, while the calculation results of the measuring points B1, C1, and D1 in the cooling process have a certain deviation compared with the test results.Figure 6
Comparison of results of the volumetric water content with numerical calculation in case 2.Figure7 shows the distribution of saturation degree in the soil for case 2 when t is equal to 1 d (heating process) and when t is equal to 5.5 d (cooling process). It can be seen from Figure 7 that the bottom saturation degree of the soil, due to gravity effect, is higher than the top saturation degree; that is, the pore water pressure at the bottom is greater. The saturation degree near the inner heat source is almost lower than that of other locations, indicating that pore water migrates from the inside to the outside due to the thermal driving effect. During the heating process, the saturation degree of the soil showed a cone-shaped distribution with a narrow top and a wide bottom, while during the cooling process, the saturation degree of the soil showed a funnel-shaped distribution with a wide top and a narrow bottom. Besides, and the direction of the pore water pressure gradient during the cooling process is from the bottom close to the boundary to the top close to the heat source.Figure 7
Saturation distribution in the soil: (a)t = 1 d during the heating process; (b) t = 5.5 d during the cooling process.
(a)(b)


## 4.1. Physical Model and Boundary Conditions

This paper uses “COMSOL Multiphysics 5.6” finite element software to solve the governing equations in Section2 and to simulate the experimental process. The calculation model adopts an axisymmetric model. The inner heat source and the outer cylindrical wall were set as impervious and incapable of displacement boundary conditions; the top cover and bottom plate were set as thermally insulated with impervious boundary conditions. The temperature boundary of the inner heat source was set to 80°C and 25°C during the rising and falling process. The temperature boundary of the outer cylindrical wall was set to 25°C, and the initial temperature of the soil was 25°C. Thus, the optimization results of relevant parameters required in the calculation are shown in Table 3. The model mesh was divided into a free triangle mesh, the maximum element size was 2.65 cm, and the minimum element size was 0.15 mm (Figure 2). The complete grid contains a total of 884 domain elements and 74 boundary elements.Table 3
Model parameters.
Parameterρs0Ea1a2ksKcpsβTβMg/cm3MPaW/(m·K)m2J/(kg·K)1/KValue1.518.321.220.681.0 × 10−229702 × 10−9−1×10−3Figure 2
The meshing of the computational model.In the calculation process, both the internal heat source temperature boundaryTin and the outer temperature boundary Tout were set as the first type of boundary conditions. As a result, the internal heat source temperature boundary is within the first 5.15 d, and the temperature was maintained at about 80°C. During the test, the temperature was basically maintained at about 25°C. It should be noted that, during the 5.15 d cooling process, the water bath needs to be replaced with cold water to achieve the cooling effect. The constant temperature water bath had no cooling function, so the internal temperature of the heat source was measured. After the temperature drop starts, it first rises and then gradually decreases to 25°C, as shown in Figure 3(a), while the outer temperature boundary has been maintained at about 25°C (Figure 3(b)), which was arranged on the outer cylindrical wall of the test device (Figure 3(b)). The temperature measured by sensors changes over time at four different locations.Figure 3
Variation of (a) heat source temperature and (b) boundary temperature with time.
(a)(b)


## 4.2. Discussion

According to the test results, it can be found that the volumetric water content of each monitoring position in the soil has changed after the temperature rising and falling process. Due to different soil-water characteristics of unsaturated soils in the process of moisture absorption and dedrying, the VG model parameters are different during the heating and cooling processes. After trial calculations, the parameterα1 = 1.89 m−1 during the heating process; it is taken as α2 = 1.93 m−1 during the cooling process. In the calculation process, the initial pore water pressure as well as pore gas pressure was 0.94 × 105 Pa and 1.01 × 105 Pa, respectively. Figure 3 shows the comparison between the calculated results and the experimental results of the volumetric water content at different monitoring positions in the soil for case 1. By comparison, it can be observed that the calculation results are consistent with the experimental results. During the heating process, the calculation results of positions B1, C1, and D1 are in good agreement with the test results, while the calculation results of position A1 are slightly lower than the test results, but the trend of change is still consistent. It can be observed that when the parameter α does not change during the heating and cooling process (α1 = α2 = 1.89 m−1) by comparing Figures 4(a) and 4(b), the calculation results present that the content of volumetric water at each location after the heating and cooling process has little change from the initial state and is far from the test result (Figure 4(a)). When the parameter α is set differently during the heating and cooling process (e.g., α1 = 1.89 m−1 and α2 = 1.93 m−1), the calculation results show that the volumetric water content at each location after the heating and cooling process is lower than the initial state and closer to the test results (Figure 4(b)).Figure 4
Comparison of the volumetric water content with numerical calculation in case 1: (a)α1 = α2 = 1.89 m−1; (b) α1 = 1.89 m−1 and α2 = 1.93 m−1.
(a)(b)Figure5 shows the comparison between the test results and the calculated results of the temperature at different monitoring positions in the soil for cases 1 and 2. In the calculation process, the parameter α of case 1 is set to α1 = 1.89 m−1 and α2 = 1.93 m−1, and the parameter α of case 2 is set to α1 = 1.89 m−1 and α2 = 1.94 m−1. By comparison, it can be seen that the numerical calculation results of the temperature changes at each measuring point during the heating process of the two cases are in good agreement with the test results, but there is a certain gap during the cooling process. This is due to the errors in the parameters and boundary conditions of the laboratory experiments and numerical calculations.Figure 5
Comparison of results of temperature with numerical calculations in case 1: (a)α1 = 1.89 m−1 and α2 = 1.93 m−1; (b) α1 = 1.89 m−1 and α2 = 1.94 m−1.
(a)(b)Figure6 shows the comparison between the test results and the calculated results of the volumetric water content at different monitoring positions in the soil for case 2. Through comparison, it can be observed that the calculation results are consistent with the experimental results. In the calculation process, the parameter α is set to α1 = 1.89 m−1 and α2 = 1.94 m−1, and the initial pore water pressure and pore gas pressure are set to 1.56 × 105 Pa and 1.46 × 105 Pa, respectively. The permeability K in saturated soil is set to 4 × 10−19 m2. It can be observed from Figure 6 that the numerical calculation results in the heating process are in good agreement with the test results, while the calculation results of the measuring points B1, C1, and D1 in the cooling process have a certain deviation compared with the test results.Figure 6
Comparison of results of the volumetric water content with numerical calculation in case 2.Figure7 shows the distribution of saturation degree in the soil for case 2 when t is equal to 1 d (heating process) and when t is equal to 5.5 d (cooling process). It can be seen from Figure 7 that the bottom saturation degree of the soil, due to gravity effect, is higher than the top saturation degree; that is, the pore water pressure at the bottom is greater. The saturation degree near the inner heat source is almost lower than that of other locations, indicating that pore water migrates from the inside to the outside due to the thermal driving effect. During the heating process, the saturation degree of the soil showed a cone-shaped distribution with a narrow top and a wide bottom, while during the cooling process, the saturation degree of the soil showed a funnel-shaped distribution with a wide top and a narrow bottom. Besides, and the direction of the pore water pressure gradient during the cooling process is from the bottom close to the boundary to the top close to the heat source.Figure 7
Saturation distribution in the soil: (a)t = 1 d during the heating process; (b) t = 5.5 d during the cooling process.
(a)(b)


## 5. Conclusion

For unsaturated porous media, the governing equations of the three-field coupling effect of thermo-hydromechanical were established based on the theory of thermodynamics and thermoelasticity. The corresponding state equations were given, and each state variable was adopted according to the temperature-dependent nonlinear expressions.A self-developed test device was utilized to discuss the hydrothermal coupling effects of two cases of unsaturated kaolin under the same temperature load. The test results were compared with the numerical calculation results. The calculation results at different positions are in good agreement with the test results, which shows that the established theoretical model is able to simulate the water-heat coupling process of unsaturated soils.Different initial saturation degrees have little effect on the change of temperature over time. Compared with the process of moisture migration, the temperature reaches a steady state faster. The characteristics of moisture migration during the temperature rise and fall process are more complicated. When the initial saturation degree of the soil is high, the volumetric water content of the measuring point during the heating process increases, and then it is in a stable state; when the initial saturation is low, the volumetric water content near the heat source increases firstly and then decreases. The volumetric water content of the soil undergoes an irreversible process during the temperature rising and falling process. The test results indicate that the volumetric water content of the measuring point after the temperature rising and falling process is lower than the initial state.



---

*Source: 2867985-2022-02-04.xml*