Abstract: This study presents a formulation for coupled moisture and air in expansive unsaturated clays. An extended critical state elasto-plastic constitutive model based on net mean stress and suction is adopted. Numerical techniques are proposed which allow for the development of an elasto-viscoplastic solution algorithm within the context of a transient finite element analysis. The representation of the development of irreversible strains within an expansive soil has been investigated via a double structure elasto-plastic relationship. The ability of such approach to simulate the accumulation of plastic strains during wetting-drying cycles has been shown.
INTRODUCTION
The engineering behaviour of unsaturated soil has been the subject of numerous experimental and theoretical investigations (Fredlund and Morgenstern, 1977; Alonso et al., 1990; Cui and Delage, 1996; Cui et al., 2002; Fleureau et al., 2002). A number of constitutive relationships have been proposed to describe the stress-strain behaviour of the material. While it is recognized that limited experimental data is available to fully characterize these relationships, much progress can be seen to have been reported (Fredlund and Rahardjo, 1993; Alonso et al., 1994, 1999; Volckaert et al., 1996; Gallipoli et al., 2003; Wheeler et al., 2003).
Suction controlled experiments (Josa et al., 1987; Alonso et al., 1995a, b; Cuisinier and Masrouri, 2004; Alonso et al., 2005) have revealed that wetting under a given confining stress may induce an irreversible volumetric deformation in unsaturated soils. An elasto-plastic approach is therefore required to represent this behaviour.
An elasto-plastic work hardening constitutive model for partially saturated soils was first proposed by Alonso et al. (1990). The deformation was related to both suction and net mean stress changes, within the context of a critical state type framework. Wheeler and Sivakumar (1995) reported results from an experimental investigation of compacted Kaolin which have provided support for the main features of the model of Alonso et al. (1990) and they suggested a number of refinements. Wheeler et al. (2003) have presented a new elasto-plastic model for unsaturated soils by including the influence of degree of saturation on the stress-deformation behaviour.
The elasto-plastic models of Alonso et al. (1990) and Wheeler and Sivakumar (1995) were developed for non-expansive and moderately expansive clays and were thought to be inappropriate for unsaturated highly expansive clays. Chu and Mou (1973) reported test results from wetting and drying cycles performed on highly expansive clay which showed a large irreversible component of swelling on first wetting. Subsequently, Alonso et al. (1995b) reported cyclic wetting-drying test results for another highly expansive clay, which showed an irreversible component of shrinkage during each drying stage. Neither of these patterns of behaviour could be represented by the existing elasto-plastic models (which assume that both wetting- induced swelling and drying-induced shrinkage are elastic processes). In the light of these experimental observations, Gens and Alonso (1992) and Alonso et al. (1994) proposed an extended form of elasto-plastic model for highly expansive unsaturated clays. The model was based on an assumption of two levels of soil fabric (an unsaturated macrofabric and a saturated microfabric within individual clay packets) with coupling between the mechanical behaviour of the two levels.
Thomas and He (1998) have presented a numerical solution using the model of Alonso et al. (1990) in which the elasto-plastic response was possible due to either conventional loading or changes in the suction in the soil. Two yield functions was then required in their formulation.
The approach presented in this study is the representation of the development of irreversible strains within an expansive soil via a double structure elasto-plastic relationship (Gens and Alonso, 1992). In this formulation, interactions between the micro- and macrostructure of swelling clays are addressed. To illustrate this approach, results of isothermal simulations using this elasto-plastic model are presented.
MATERIALS AND METHODS
Theory and Formulation
Unsaturated soil is a three-phase porous medium consisting of solid,
liquid and air. The deformation of unsaturated soil is coupled to both
stress variation and moisture movement. An extended critical state model
proposed by Alonso et al. (1994) is adopted here to relate the
deformation in the variations of stress and suction in the soil.
Stress-Strain Constitutive Relationship
A double structure elasto-plastic model based on the framework for
expansive soils presented by Gens and Alonso (1992) is described. This
framework assumes that a double structure model can describe the behaviour
of expansive soils. The double structure is defined as consisting of a
macro- and a microstructure. A slightly swelling model (Alonso et al.,
1990) is utilised to describe the macrostructure of the soil, defined
as the level at which the major structural arrangements take place. The
section of the model describing the microstructure, defined as the level
at which basic swelling of active minerals occurs, assumes a saturated
state. The authors propose that an effective stress law may be applied.
Only reversible strains are allowed to develop in the microstructure whilst
irreversible strains can develop in the macro level. Coupling of the two
levels allows only the microstructure to influence the macrostructure.
This coupling provides a mechanism that allows the development of strains
in the microstructure to lead to softening/hardening of the macrostructure
by the development of irreversible strains. This is accomplished by assuming
that the irreversible strains will occur when the effective stress goes
beyond either a maximum or a minimum level. The authors propose that these
maximum and minimum levels of effective stress be defined by a pair of
yield surfaces between which the stress state of the soil lies. These
are known as the suction increase SI and the suction decrease SD yield
surfaces, respectively. Furthermore, it is assumed that these two lines
are coupled, leading to the movement of either when a load path move to
the other.
Following the approach presented by Gens and Alonso (1992), the increment of total strain can be defined as:
(1) |
where, the subscripts M and m represent the macro and micro contributions, respectively and the superscripts e and p represent elastic and plastic contributions.
The expression for the macrostructural elastic volumetric strain component due to stress changes is given by Alonso et al. (1990) as:
(2) |
is adopted here. κ is the stiffness parameter for changes in net mean stress p in elastic region and v is the specific volume. Alonso et al. (1990) also defined the macrostructural elastic volumetric strain component due to suction changes, at a constant stress, as:
(3) |
where, κs is the stiffness parameter for changes in suction in elastic region, s is suction and patm is the atmospheric pressure.
As discussed earlier, the framework is based on the assumption that the microstructure remains saturated even when the macrostructure has a negative pore water pressure and that the principle of effective stress will therefore be valid within the aggregates. In term of net mean stress and suction, effective mean stress is defined as (p+s). The microstructural elastic volumetric strain component dεem is therefore expressed as a function of (p+s), that is (Alonso et al., 1994):
(4) |
where, αm and βm are material parameters and A(p+s) is the microstructural stiffness.
The yield function related to the preconsolidation stress, for isotropic stress states is identical to the yield function presented by Alonso et al. (1990) for the slightly swelling model, i.e.,
(5) |
where, p is the net mean stress, q is the deviatoric stress and p0 is the preconsolidation stress, ps is the parameter controlling suction effect on cohesion and M is the slope of the critical state line.
The preconsolidation stress p0, was defined as:
(6) |
(7) |
where, p*0 is the preconsolidation stress of saturated soil, pc is the reference stress, λ(0), λ(s), β and r are the parameters controlling the stiffness of the soil.
In addition to the Load Collapse (LC) yield surface, two further equations are required to define SI and SD yield surfaces. Based on the above assumptions of the definition of these surfaces, Alonso et al. (1994) presented them as:
(8) |
(9) |
where, sI and sD are stress parameters that can physically be defined as the level of suction at which, under zero mean confining stress, irreversible shrinkage or swelling strains are produced on drying or wetting paths. The authors also present two yield loci to define the amount of plastic strains produced when the surfaces yield. From these, two hardening laws for the SI and SD yield surfaces were presented by Alonso et al. (1994) as:
(10) |
(11) |
where, hI and hD are functions defining the relationship between microstructural elastic and plastic strains. With suitable substitution of Eq. 1 the stress-strain relationship can be expressed as:
(12) |
which can be rewritten as:
(13) |
where, σ” is the net stress.
The above constitutive relationship has been incorporated in a more general numerical model to achieve a solution satisfying the stress equilibrium equation:
(14) |
where, δij is the Kronecker delta function and bi is the body force.
Moisture Transfer
The velocity of liquid is based on the generalised Darcy’s law and
is then expressed as:
(15) |
where, Vw is the velocity of the liquid, Kw is the non saturated hydraulic conductivity, μw is the pore water pressure, γw is the unit weight of water et z is the reference position and ∇ is the gradient operator.
An expression of the degree of saturation is obtained from the experimental results performed on Boom clay by Centro International de Metodos Numericos en Ingeniera, UPC, Spain, which gives:
| |
(16) |
where, s is the suction in Pa.
The hydraulic conductivity relationship for Boom clay is expressed by Alonso
et al. (1995b) as:
(17) |
The governing equation for moisture transfer in an unsaturated soil can be expressed as:
(18) |
where, ∇ is the divergence operator
Dry Air Transfer
Air in unsaturated soil is considered to exist in two forms: bulk
air and dissolved air (Fredlund and Rahardjo, 1993). In this approach
the proportion of dry air contained in the pore liquid is defined using
Henry’s law.
The pore air velocity is assumed to satisfy the Darcy’s law, i.e.,
(19) |
where, Ka is the air conductivity, ua is the pore air pressure and γa is the unit weight of air.
Alonso et al. (1995b) have proposed a relationship for air conductivity for Boom clay as:
(20) |
where, μa is the air viscosity.
To take into account the dissolved air in water, Henry’s law is used and the governing equation of the air flow can be expressed as (Thomas and Sansom, 1995):
(21) |
where, Hs is Henry’s volumetric coefficient of solubility and ρa is the density of air.
Numerical Algorithm
The complexity and coupled nature of the governing differential equations presented in previous sections prevents a direct analytical solution. A numerical approach is then required and the formulation of this solution is based on the finite element method for the spatial discretisation and on a backward difference mid-interval time stepping scheme for temporal discretisation. Before the elasto-plastic relationship can be developed, the plastic contribution needs to be addressed. A viscoplastic algorithm has been found to give efficient solutions to elasto-plastic problems (Thomas and Harb, 1990).
In this study the Galerkin weighted residual approach (Zienkiewicz and Taylor, 2000) is employed. Two dimensional eight noded isoparametric elements are used. The system of the three Eq. 14, 18 and 21 can be discretised spatially and converted as:
(22) |
uw, ua and u are pore water pressure, pore air pressure and displacement, respectively. The dot above variables refers to the derivative with respect to time. Kij et Cij represent the corresponding matrices of the governing equations with (I, j = w, a, u).
For simplicity, it is convenient to rewrite Eq. 22 in the following form:
(23) |
where, {φ} resent the global unknown vector, i.e., {uw ua u}T.
To solve Eq. 23, a general form of fully implicit mid-interval backward difference time stepping algorithm is used to discretise the governing equation temporally. Therefore Eq. 23 can be rewritten as:
(24) |
(φn) is the level of time at which the matrices K, C and f are to be evaluated and is given by:
(25) |
where, ω is integration factor, which defines the required time interval (ω ε [0, 1]) and θ = 0, 0.5, 1 for backward, central and forward difference schemes, respectively.
For a mid-interval backward difference scheme, ω = 0.5 and θ = 0. Therefore, Eq. 24 reduces to:
(26) |
Equation 26 may be rewritten in alternate form as:
(27) |
A solution for {φn+1} can be obtained provided the matrices K, C and f at time interval n+0.5 can be determined. This is achieved by the use of a predictor-corrector iterative solution procedure. The algorithm of this procedure is presented below:
| • | Evaluate K, C and f at time n to produce a first estimate named the predictor. |
| • | Use the predictor and previous time step values to evaluate K, C and f at time n+0.5, producing an estimate called corrector. |
| • | Check for yielding, if so, evaluate the viscoplastic
strain rate |
| • | Check for convergence of successive correctors via either of the following conditions: |
(28) |
(29) |
where, i is the iterative level, c represents the case of corrector value and TLabs and TLrel are the matrices of absolute and relative tolerances for each variable. Also, check stress equilibrium to ensure the residual force is within a tolerance value. The residual force Ψ can be calculated (Owen and Hinton, 1980) as:
(30) |
where, ΔF is the increment of applied force and B is the strain displacement matrix. If both of the variable convergence conditions are not met or stress equilibrium is not satisfied then return to point 2. At this point the corrector becomes the predictor.
| • | Move onto next time step and repeat from point 1. |
The number of corrector iterations required for convergence is related to the size of the time step used. This property has been used to modify the time step size used as a simulation processes. If the number of iterations exceeds a specified value, the time step size is reduced by a factor. Likewise, if the number of iterations required falls below a specific value, the time step is then increased by a factor. This variable time stepping scheme helps to provide an efficient solution algorithm.
RESULTS AND DISCUSSION
Numerical Example
The example presented in this study is a simulation of a set of wetting-drying
cycles in a suction controlled oedometer tests on Boom clay performed
by Alonso et al. (1995a). First, the soil sample was subjected
to an increase in vertical stress to the required level without lateral
deformation. Then it was soaked to saturation while the vertical stress
was kept constant. The test was performed under various net vertical tress
levels. It was observed that the soil sample with lower initial net stress
produce considerable swelling under wetting condition. However those soil
samples that had initial net stress higher than 200 kPa showed significant
wetting collapse behaviour after an early swelling stage during the test.
The higher the net stress level, the more wetting collapse deformation
is observed. When the net vertical stress was increased to over 300 kPa,
the soil sample did not show any swelling during the whole wetting path.
The experimental results provide therefore a good example which displays
the typical feature of unsaturated soil.
The domain geometry of size 10 mm by 100 mm was discretised into 10 eight noded isoparametric elements equally sized (10 mm by 10 mm). The results for this mesh have been found to give spatially converged results when compared to a finer mesh. The deformation was only allowed to develop in the vertical direction. All the parameters and functions used in the formulations are shown in Table 1.
Figure 1-5 show the deformation response for the five tests reported. With suitable coupling between SI and SD yield surfaces it has been possible to simulate the accumulation of plastic strains during the wetting-drying cycles.
| Table 1: | Material parameters used in the numerical example |
| Fig. 1: | Volumetric strain due to change of suction under net vertical load of 19.62 kPa |
| Fig. 2: | Volumetric strain due to change of suction under net vertical load of 49.05 kPa |
| Fig. 3: | Volumetric strain due to change of suction under net vertical load of 98.01 kPa |
| Fig. 4: | Volumetric strain due to change of suction under net vertical load of 206.01 kPa |
| Fig. 5: | Volumetric strain due to change of suction under net vertical load of 392.40 kPa |
The results show the ability of a double structure approach to qualitatively represent the trends of behaviour observed in the samples.
CONCLUSION
The research presented in this study describes recent research performed on the development of a model to describe the coupled flow of moisture and air in expansive soils. The representation of the development of irreversible strains within an expansive clays has been investigated via a double structure elasto-plastic relationship. The simulation performed has produced a set of very encouraging results. The model has been able to represent the experimentally measured trends of behaviour of the sample during initial wetting and drying path and the subsequent scanning cycles.