Finite element analyses of rate-dependent thermo-hydro-mechanical behaviors of clayey soils based on thermodynamics

Abstract

Clayey soils in the vicinity of energy geostructures may be exposed to long-term periodic thermal cycles. The creep and consolidation behaviors of the clayey soils can be both rate-dependent and temperature-dependent, and the underlying physical mechanisms are merely investigated theoretically. In this study, based on the theory of thermodynamics, a fully coupled thermo-hydro-mechanical (THM) finite element (FE) program for saturated soils is developed for this purpose. The FE formulation accounts for the combined effect of rate and temperature through the novel concept of granular temperature. Simulations of THM coupled validation cases and a series of experimental observations on the soft Bangkok clay are carried out. The obtained numerical results exhibit good agreement with analytical solutions and laboratory measurements. It is found that three fundamental physical mechanisms contribute to the irreversible thermal contraction observed for normally consolidated and lightly overconsolidated clays under drained thermal cycles: (1) the thermal creep excited by mass exchange from adsorbed water to free water; (2) the mechanical creep induced by confining stresses; and (3) the increase in granular packing caused by the thermal expansion of soil particles. The thermal contraction generally stabilizes within a few thermal cycles, as a result of the noticeable reduction in the thermal creep rate. It is further demonstrated that the transient heat transfer and the heating rate can greatly influence the deformation of clays subjected to thermal cycles.

Appendices

Appendix 1: Detailed expressions of the matrices for the THM coupled FE formulation

$$ \left[ {D_{\varepsilon } } \right] = \left[ {\begin{array}{*{20}c} {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \varepsilon_{r}^{e} }}} & {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \varepsilon_{\theta }^{e} }}} & {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \varepsilon_{z}^{e} }}} & {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \varepsilon_{rz}^{e} }}} \\ {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \varepsilon_{r}^{e} }}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \varepsilon_{\theta }^{e} }}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \varepsilon_{z}^{e} }}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \varepsilon_{rz}^{e} }}} \\ {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \varepsilon_{r}^{e} }}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \varepsilon_{\theta }^{e} }}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \varepsilon_{z}^{e} }}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \varepsilon_{rz}^{e} }}} \\ {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \varepsilon_{r}^{e} }}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \varepsilon_{\theta }^{e} }}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \varepsilon_{z}^{e} }}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \varepsilon_{rz}^{e} }}} \\ \end{array} } \right] $$

$$ \left[ {D_{T} } \right] = \left[ {\begin{array}{*{20}c} {\frac{{\partial \sigma_{r} ^{\prime}}}{\partial T}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{\partial T}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{\partial T}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{\partial T}} \\ \end{array} } \right]^{T} $$

$$ \left[ {D_{\rho } } \right] = \rho_{d} \left[ {\begin{array}{*{20}c} {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{r} ^{\prime}}}{{\partial \rho_{d} }}} & 0 \\ {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{\theta } ^{\prime}}}{{\partial \rho_{d} }}} & 0 \\ {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \sigma_{z} ^{\prime}}}{{\partial \rho_{d} }}} & 0 \\ {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \rho_{d} }}} & {\frac{{\partial \tau_{rz} ^{\prime}}}{{\partial \rho_{d} }}} & 0 \\ \end{array} } \right] $$

$$ \left\{ {\dot{Y}} \right\} = \left[ {\begin{array}{*{20}c} {(T_{g} )^{0.5} \varepsilon_{r}^{e} + (m_{1} - \frac{1}{3})(T_{g} )^{0.5} \varepsilon_{v}^{e} } \\ {(T_{g} )^{0.5} \varepsilon_{\theta }^{e} + (m_{1} - \frac{1}{3})(T_{g} )^{0.5} \varepsilon_{v}^{e} } \\ {(T_{g} )^{0.5} \varepsilon_{z}^{e} + (m_{1} - \frac{1}{3})(T_{g} )^{0.5} \varepsilon_{v}^{e} } \\ {(T_{g} )^{0.5} \gamma_{rz}^{e} } \\ \end{array} } \right] $$

$$ [K] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {B^{T} ([D_{\varepsilon } ] + [D_{\rho } ])Brdrdz \, }} $$

$$ \left[ {K_{p} } \right] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {B^{T} \left\{ m \right\}N_{p} rdrdz \, }} $$

$$ \left[ {K_{T} } \right] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {B^{T} [D_{T} ]N_{T} rdrdz \, }} $$

$$ \left\{ F \right\} = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( { - B^{T} [D_{\varepsilon } ]\left\{ {\dot{Y}} \right\} + N^{T} \left\{ {\dot{g}} \right\}[\rho_{d} + \rho_{l} (1 - \frac{{\rho_{d} }}{{\rho_{s} }})]} \right)rdrdz}} + \sum\limits_{{S_{\sigma }^{e} }} {\int\limits_{{S_{\sigma }^{e} }} {N^{T} \left\{ {\dot{s}} \right\}rdS} } $$

$$ [H_{\varepsilon } ] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( {A_{1} N_{p}^{T} m^{T} B} \right)rdrdz \, }} $$

$$ [H] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( {A_{2} N_{p}^{T} N_{T} } \right)rdrdz \, }} $$

$$ [H_{p} ] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( {\frac{\kappa }{\mu }A_{1} (\nabla N_{p} )^{T} \nabla N_{p} } \right)rdrdz \, }} $$

$$ [H_{T} ] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left[ {\frac{\kappa }{\mu }\beta_{l} \left( {N_{p}^{T} \left\{ p \right\}^{T} (\nabla N_{p} )^{T} \nabla N_{T} - \rho_{l} N_{p}^{T} \left\{ g \right\}^{T} \nabla N_{T} } \right)} \right]rdrdz \, }} $$

$$ \left\{ {F_{p} } \right\} = \sum\limits_{{S_{p}^{e} }} {\int\limits_{{S_{p}^{e} }} {\left[ {A_{1} N_{p}^{T} \cdot (\tilde{v}_{r} + \tilde{v}_{z} )} \right] \, rdS} } + \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left[ {\frac{\kappa }{\mu }A_{1} \rho_{l} (\nabla N_{p} )^{T} \left\{ g \right\}} \right]rdrdz \, }} $$

$$ [R] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( {\sum {\rho_{\alpha } \varphi_{\alpha } c_{\alpha } } } \right)N_{T}^{T} N_{T} rdrdz \, }} $$

$$ [R_{T} ] = \sum\limits_{{\Delta^{e} }} {\iint\limits_{{\Delta^{e} }} {\left( {\lambda \left( {\nabla N_{T} } \right)^{T} \nabla N_{T} + \rho_{l} \varphi_{l} c_{l} \frac{\kappa }{\mu }\left[ {N_{T}^{T} \left\{ p \right\}^{T} \left( {\nabla N_{p} } \right)^{T} \nabla N_{T} - \rho_{l} N_{T}^{T} \left\{ g \right\}^{T} \nabla N_{T} } \right]} \right)rdrdz \, }} $$

$$ \left\{ {F_{T} } \right\} = - \sum\limits_{{S_{T}^{e} }} {\int\limits_{{S_{T}^{e} }} {\left[ {N_{T}^{T} \cdot (\tilde{\eta }_{r} + \tilde{\eta }_{z} )} \right]rdS} } $$

Appendix 2: THM processes in saturated linear elastic geomaterial

The one-dimensional THM processes in this case study was first proposed by Suvorov and Selvadurai [52] for analyzing the quasi-static response of saturated linear elastic geomaterials subjected to heating. The case serves as a valuable validation of numerical solutions because a rigorous analytical solution has been provided [52].

Consider a one-dimensional elastic soil column of 10 m in height, as shown in Fig. 

Fig. 12

FE mesh and boundary conditions of the THM processes in the saturated elastic soil column

12. At time t = 0, the column is instantaneously heated from 0 to 100 °C without drainage. During the time ,the bottom end is fixed while the top end has no displacement limit. Accordingly, excess pore pressure builds up within the soil column and the top end elevates due to thermal expansion. Once the entire column reaches 100 °C, the pore pressure and the temperature at the top end are immediately reduced to zero and a constant compressive stress of 10 MPa is applied to the top. Temperature, pore pressure and displacement are observed within the next 365 days.

Figure 12 shows the geometry, the FE mesh and the boundary conditions of the problem. Table

Table 2 Soil properties for the validation case in "Appendix 2"

2 lists the soil properties used in the analytical solution by Suvorov and Selvadurai [52] and the way to convert the properties into the FE code. The heat convection is not considered in this case. Degradation from the extended GSH model to a linear elastic model can be completed using the following equation, provided that all other hyperelastic and plastic constitutive parameters are set to zero

$$ \left\{ \begin{gathered} B_{0} = \frac{5}{12}\frac{E}{(1 - 2\nu )} \hfill \\ \xi = \frac{E}{{2B_{0} (1 + \nu )}} \hfill \\ \end{gathered} \right. $$

(20)

where E and v are the Young’s modulus and the Poisson’s ratio.

Figure 

Fig. 13

Comparison between the numerical and the analytical solution for the evolution of a temperature, b pore pressure and c displacement

13 shows the evolution of temperature, pore pressure and displacement with time and depth in 365 days. A perfect match is found between the FE calculation and the analytical solution. Since the heat convection is not considered, the temperature evolution in Fig. 13a is totally controlled by the effective specific heat and the thermal conductivity according to Eq. (17). In Fig. 13b, the initial excess pore pressure of 36.3 MPa is induced by the elevated temperature (which contributes to 26.3 MPa) and the 10 MPa compressive stress. At Day 100, not only the positive pore pressure induced by heating has been dissipated, but a negative pore pressure has been generated. This is because the thermal expansion coefficient of the liquid phase is set to zero while that of the solid phase isn’t. During the cooling stage, the contraction of solid phase adds a tensile force to the liquid phase which results in the negative pore pressure. Figure 13c shows the displacement along the elastic soil column. The initial displacement is caused by the thermal expansion. The final displacement after temperature returns to zero depends on the relationship between the soil stiffness and the compressive stress applied. In general, the simulation proved the accuracy of the FE code in dealing with THM coupled elastic problems.