Geofluids

Volume 2018, Article ID 4089612, 14 pages

https://doi.org/10.1155/2018/4089612

## Numerical Analysis Method considering Coupled Effects of THMC Multifields on Unsaturated Expansive Soil Subgrade Treated with Lime

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China

Correspondence should be addressed to Zheng Lu; moc.361@msrhwzl

Received 29 June 2017; Revised 29 January 2018; Accepted 11 February 2018; Published 21 March 2018

Academic Editor: Qinghui Jiang

Copyright © 2018 Jie Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The response model subjected to coupled effect of thermo-hydro-mechano-chemical (THMC) was built in the context of basic theories in the polyporous polyphasic medium mechanics, the mixture theory in the continuum mechanics, and the thermodynamic theories. The finite element discretization of the response model was implemented based on the Galerkin method. The processes of salt leaching and accumulating were analyzed in the numerical results. The results among the numerical results and measured results were compared and discussed. Finally, the solutes migration rule of the soil subjected to the atmosphere eluviations was revealed, and the reasonableness of the coupling model and the finite element method was proved. The agreement between the numerical and the measured results was good, which indicates that the THMC model and finite element program were useful in solving the coupling problems of unsaturated soil. Moreover, the salt dissolution process has a larger effect on the salt movement compared to that of the salt accumulation process. Comparing with the salt leaching effect caused by rainfall, the salt accumulation effect caused by evaporation was smaller.

#### 1. Introduction

In China, the stability of subgrade is the basic premise to guarantee highways becoming the leading infrastructure construction. Expansive soil cannot be directly used as the subgrade filling within the subgrade. Engineering properties and mechanical strength indexes of the expansive soil can be increased significantly through modification using the lime. However, the expansive soil subgrade treated by the lime will easily suffer from repetitive drying-wetting cycles and unceasing strength degradation in the rainy and hot weather, which would make the subgrade failure. It is essential to study the effects of solute migration in the expansive soil treated by lime under the atmospheric actions on the subgrade stability. A lot of domestic and overseas experts have developed related research about soil subgrade treatments by lime. Kumar et al. [1] studied, compared, and analyzed variations of compaction and strength properties of the expansive soil before and after treatment by the lime. Nowamooz and Masrouri [2] researched volumetric strains caused by variations in suction or stress of the expansive mixture of limestone soil and flyash. Guney et al. [3] studied the expansibility of expansive soil treated by the lime under drying and wetting cycles. Estabragh et al. [4] researched mechanical behaviors of the expansive soil experienced drying and wetting cycles with different wetting fluids. Tripathy and Subba Rao [5] studied the cyclic swelling and shrinking behaviors of tamped expansive soil through theoretical analysis and indoor experiments. Rajeev et al. [6] researched variation trends of humidity and temperature of the subgrade in the long term under atmospheric actions. Al-Mukhtar et al. [7, 8] studied the micromechanism and mechanical behaviors of the expansive soil treated by the lime. Wang et al. [9] researched variation characteristics of kinetic parameters of the expansive soil treated by the lime. Stoltz et al. [10] performed multiscale analysis of swelling and shrinking of the expansive soil treated by the lime from the microscopic view. Ojuri et al. [11] studied the geotechnical and environmental evaluation of lime-cement stabilized soil-mine tailing mixtures for highway construction. Dash and Hussain [12] researched the influence of lime on shrinkage behavior of soils. Hotineanu et al. [13] studied the effect of freeze-thaw cycling on the mechanical properties of lime-stabilized expansive clays.

The above researches show that a lot of experts have performed fruitful researches on physical and mechanical properties of the expansive soil treated by the lime from macroscopic and microscopic views. However, considering atmospheric actions, soil modification, and solute migration rule to study the long-term properties of the expansive soil treated by the lime in the highway subgrade would got closer to the practical results. In this paper, as for the expansive soil treated by the lime under atmospheric actions, the thermo-hydro-mechano-chemical coupled model was built on the basis of the Polyporous Polyphasic Medium Mechanics, the mixture theory in the continuum mechanics, and thermodynamic theories of the irreversible process. The finite element discretization of the response model is implemented based on the Galerkin method. The processes of salt leaching and accumulating are analyzed in the numerical results. The numerical results and measured results are compared and discussed. Finally, the solutes migration rule of the soil subjected to the atmosphere eluviation is revealed, and the reasonableness of the coupling model and the finite element method is proved.

#### 2. Coupled Equations for THMC

##### 2.1. Basic Assumptions

Some basic assumptions were given before building THMC model to make the problem has universality: (a) the unsaturated soil is polyporous medium, its basic framework is particles, and liquid water and gas exist in the pores of the particle framework. (b) The gas phase is ideal gases with the mixture of dry air and water vapor. The liquid phase consists of water and gases dissolving in water. (c) Mechanical behavior of the polyporous medium (gases, liquids, and solids) satisfies the principle of local stress balance, and the entire framework satisfied the principle of effective stress. (d) The temperature of the soil is the average temperature, which satisfied the local equilibrium hypothesis, heat transfer satisfies the Fourier law, and the framework and water can be considered as isotropic thermoelastic materials. (e) The unsaturated soil satisfies general assumptions and the continuum hypothesis of the mixtures theory. (f) Actions of degradation and stagnant water effect are not taken into account in the solute equilibrium equation.

##### 2.2. Balance Equations

In the saturated-unsaturated soil mass, the total stress increment and effective stress increment can be expressed in the incremental forms as follows, respectively:where is the tangent elastoplastic modulus matrix; is the total strain increment; is the strain increment caused by the pore pressure and ; is the strain increment caused by temperature, and ; is the strain increment caused by swelling, and ; is the unit array of normal stress; is the thermal expansion coefficient of solid particles; is the water swelling expansion coefficient; is the bulk modulus of solid particles; is the average stress of gases and liquids inside pores, and ; and are saturability of the pore water and pore gas, respectively, and ; and are pressures of the pore water and pore gas, respectively.

The stress equilibrium equation is stated as follows:where is the body force.

According to the principle of virtual work, the incremental stress equilibrium equation can be stated as follows:where is the total stress increment; and are increments of body force and surface force, respectively; and are virtual strain and displacement, respectively.

Equations (2) and (3) are substituted into (4), and the following is obtained:According to the law of conservation of mass, the fluid continuity equation can be obtained:Similarly, the equation of solid mass conservation can be obtained:Equation (7) is the mass conservation equation of the component in the phase. It is the liquid phase as is equal to . It is the gas phase as is equal to . The component is water as is equal to , and the component is the gas as is equal to . For instance, and indicate air resolved in the liquid phase and water in the liquid phase, respectively; and indicate dry air and water steam in the gas phase, respectively; is the decreasing or increasing rate of unit volume water during evaporation or condensation. and are absolute moving velocities of solid particles and fluid, respectively; is saturability of the water phase or the gas phase.

The fluid flux , which flows across the unit area per unit time, is stated as follows:where is the diffusion flow of component in the phase; is the relative flow when there are movements of the solid phase in the deformed porous medium, and the subscript indicates that the velocity and flow are relative to the solid phase. According to the definition of relative velocity, it can be obtained that . Equation (6) can be stated as follows:According to (6), (7), and the definition of material derivative, , it can be obtained thatSubstituting (10) into (11), the equilibrium equation of phase component in the deformed porous medium can be obtained as follows:Under the situation of small deformation, items such as and can be neglected; can be presented as the volumetric strain of the framework . Equation (12) can be simplified as follows:In (13), the right item is the material flow including convective (Darcy) and nonconvective diffusion; the left items indicate the change rates of fluid material caused by the saturability and density (the first item), volumetric strain (the second item), volumetric deformation of soil particles (the third item), and phase change (the fourth item), respectively.

Considering the mass conservation of every component, , the total equilibrium equation is as follows:As for the water including liquid water and steam in the air (), the total equilibrium equation is as follows:Equation (14) and (15) are continuity equations of the gas and liquid, respectively. For the sake of convenience in writing, in the following derivation process, , , , and .

Variations of densities of steam, liquid water, and particles with time can be expressed aswhere ; is the density of fresh water; is the density ratio of solute to water in the pore water; is the compression modulus of fresh water; is the thermal expansion coefficient of fresh water; is the compression modulus of particles; is the thermal expansion coefficient of the framework; and is the Biot coefficient.

Considering effects of temperature and solute concentration in the unsaturated soil, the effusion flux of water steam can be described using the revised Fick law as follows:where is the effusion flux of water steam, with a unit of kg/m^{2}/s; is the factor considering pore curvature in the soil [14], and is steam diffusion coefficient [15].

Considering the effects of solute in the unsaturated soil, the density of steam can be expressed as

Substituting (16)~(18) into (15), the continuous equilibrium equation of pore water considering combined effects of framework deformation, steam, solute, and temperature can be obtained as follows:The coefficients in the equation are shown in Appendix A.

##### 2.3. Energy Conservation Equation

The total internal energy of every fluid phase and energy exchange flux can be expressed asSimilarly, as for the solid phase, it can be obtained:where is a unit mass of internal energy of phase ; and are energies of water and gas per unit mass.

According to the total energy conservation and the assumption of local thermal equilibrium, it can be obtained:where is the average heat conduction flux; , , and are convection fluxes of the liquid phase, gas phase, and solid phase, respectively; is the term of source and convergence. Effects of thermal conduction in deformation of solid phases are neglected.

When variations of temperature and pressure are not large enough, the internal energy in the energy balance equation can be replaced by the specific enthalpy. Specific enthalpies of the solid phase, water, and water steam can be expressed as follows:where , , and are specific enthalpies of the solid phase, water, and water steam, respectively; and are the specific heat of solids and water; is the specific heat at constant pressure of water steam; is the latent heat of vaporization of water steam.

Replacing the internal energy with enthalpy in (23a) and (23b),

In the above derivation, effect of the density of solids on energy is neglected, and it is assumed that the specific heat can be considered as the invariant when variations of pressure and temperature are not large.

The right part of the energy equation (22) can be abbreviated as follows:

By substituting (23a), (23b), and (27) into (22), the energy conservation equation can be obtained asThe coefficients in the equation are shown in Appendix B.

##### 2.4. Solute Transport Equation

When the temperature does not vary significantly, it can be considered that migration of solute under the temperature gradient can be neglected; therefore the solute transport equation can be obtained according to the mass conservation law of the solute:where is the amount of adsorbed solute in soil per mass, which is a function of concentration. Then (30) can be expressed aswhere is the volumetric water content of pore water; is the density of dry soil; is the hydrodynamic dispersion coefficient of the solute; is the water flux; is the moisture capacity.

##### 2.5. Initial Conditions and Boundary Conditions

Initial conditions inside the given analysis area and on the boundary when is zero are as follows:

Dirichlet boundary conditions can be expressed as follows:

Neumann boundary conditions can be expressed as follows:where , , , , and are assigned liquid water flow, water steam flow, solute flow, and tensile force on the boundary . and are the mass density and temperature of water steam at undisturbed locations far away from boundaries; is the mass convection coefficient; is the thermal convection coefficient. Mixed boundary conditions may appear during the actual calculation. Therefore boundary conditions should be determined according to the actual situation.

The initial concentration distribution is given aswhere is the depth of calculated soil layer.

At the earth surface where is zero, boundary conditions of solute migration can be one kind of the Dirichlet, Neumann, and Cauchy boundary conditions.

(1) As the solute concentration at the earth surface is given, the Dirichlet boundary condition is expressed as follows:(2) As the earth surface is in the status of infiltration and the solute concentration or the degree of mineralization of rainfall or irrigation water is given, it is the Cauchy boundary condition expressed as follows since the solute transport flux :

As the intensity of water supply is smaller than the soil infiltration capacity, the water movement flux at the earth surface in the above equation is equal to the intensity of water supply ; as the intensity of water supply is larger than the soil infiltration capacity, the water movement flux at the earth surface can be calculated through the above equation.

(3) As the earth surface is in the status of evaporation, the intensity of evaporation adopts the Penman-Wilson formula of the actual evaporation amount at the surface of the unsaturated soil as the evaporation boundary condition:where is the slope of the vapor pressure-temperature curve of the saturated steam; is the amount of net radiation at soil surface; is the humidity constant; is the amount of potential evaporation and . Here is a function of wind, and , where is the wind velocity; is the vapor pressure of air above the evaporating surface; is the reciprocal of relative humidity of the air, which is ; is the reciprocal of relative humidity of the soil surface, which is .

Here it belongs to the Cauchy boundary condition, which can be expressed as follows since the solute flow is zero at the earth surface during evaporation:(4) During field or indoor experiments, as it has neither evaporation nor infiltration on the earth surface, it belongs to the Neumann boundary condition, and it can be obtained:

##### 2.6. Finite Element Discretization of the THMC Coupled Governing Equation

Analytical solutions of the THMC coupled mathematical model can be hardly achieved due to the complexity of solution conditions and irregularity of solution region. Therefore, the only way to solve the model is to obtain approximate solutions through numerical methods. The THMC coupled governing equations such as (5), (19), (29), and (31) are spatial discretized as follows:whereEquation (36) can be expressed in the form of matrix as follows:

The Monolithic Augmentation Approach is applied to solve the model. Changes within every time step are assumed as linear, which meanswhere is a point within the time step ; the item at the left of the equation contains the values of , , , and at the time of , which means ; and .

Therefore (38) can be transformed into the following iterative solution equation:

The above equation is the spatial discretization of nonlinear equations. Hence the nonlinear governing equation is transformed into linear problem through discretization of space and time. The asymmetric stiffness matrix can be directly solved by using the asymmetric solver.

#### 3. Assessment and Validation of the Numerical Model

##### 3.1. Simulation of Salt Leaching and Accumulating

Combined with the THMC coupled governing equation and initial conditions and boundary conditions (34a)~(34b) and (35a)~(35f), as well as the given initial conditions and boundary conditions of water flow and temperature, the mathematic model of solute migration problem in the soil treated by the lime in the subgrade was then built.

We choose a calculation example to study the salt leaching and accumulating processes in field to investigate the migration rule of salinity in the unsaturated soil, and the applicability and accuracy of the model are verified: dynamic test was developed to study the changes of the water and salt of the soil in the condition of fresh water and salt water flushing [16]. The test was conducted on a plot (6.1 m × 6.1 m). The average water content of the soil was 0.2. The observation and sample point was set in the depth of 30, 60, 90, 120, 150, and 180 cm, respectively. Based on the test, the calculated results of the THMC model we raised will be contrasted with the test results, and the model will be verified.

According to the observation data, Bresler [17] has calculated the experimental process above, and the soil water characteristic curve is generalized as follows:

The relationship between hydraulic conductivity and water content is

During calculation, the hydrodynamic dispersion coefficient adopts the following equation:where the value of is 0.04 cm^{2}/min; ranges from 0.001 to 0.005; equals 10; ranges from 0.28 to 0.55 and equals 0.28, 0.39, and 0.55 in the calculation, respectively.

Since there are no calculating parameters such as temperature and volumetric deformation, the model was simplified as the water-solute coupled model, and the dynamic of soil water and salt was numerically simulated using the Bresler hydraulic characteristic parameters. Figure 1 showed the comparison between numerically calculated results and measured results, which indicated that distributions of water content and concentration were consistent with measured data. After two hours of irrigation using the saline solution, salinity migrated downward with water at a depth up to 30 cm, and the concentration decreased when the depth increased from zero (at the earth surface) to 30 cm. Then the salt was leached using the fresh water. After nine hours and eleven hours of salt leaching, wetting fronts occurred at depths of 35 cm and 110 cm, respectively, while solute fronts (the maximal concentration) appeared at depths of 40 cm and 50 cm, with peak values decreased from 210 to 170 and 150, respectively. After thirteen hours of salt leaching, solute front appeared with a depth of 60 cm with a peak value of 140. After seventeen hours of salt leaching, wetting front appeared at a depth of 170 cm.