c
TU Delft, The Netherlands, 2006
IMPLEMENTATION AND EVALUATION OF A BINARY
MIXURE MODEL AND FE-ANALYSIS ON WATER
SATURATED SOIL
Aylin Ahadi
Lund University, Division of Mechanics Box 118, S-221 00 Lund, Sweden e-mail:Aylin.Ahadi@mek.lth.se
Key words: binary mixture model, numerical simulation, FEM, water saturated, soil Abstract. A binary mixture model for soils has been formulated and analyzed. In the
model one of the constituents is soil and the other constituent an incompressible liquid. The effects of water are included by formulating a mixture model is formulated as an ex-tension of an existing non-associated elasto-plastic material model for granular materials. This paper is focused on the balance of momentum equations of the mixture model, pri-marily to avoid the difficulties arising when the interaction terms between the constituents needs to be formulated and constituted. Considering the momentum balance of the mixture as an extension of the balance equation of the soil constituent, it is supplemented with body forces, stresses and flow effects from the water constituent. These supplemented terms are evaluated and analyzed. An alternative way of deriving equation of motion for the fluid is also presented. The computer algorithm for integrating the constitutive relation of dry sand is supplemented and numerical simulations with a finite element code are carried out. It is shown that although the mixture model consists of two constituents, it is possible to include the existence of voids and a system with three constituents, i.e. water unsatu-rated soil, can be analyzed. Results from single element tests and simulations of a flexible footing resting on a surface of water saturated sand are presented.
1 INTRODUCTION
incompressible Navier-Stokes fluid. To avoid the difficulties arising from formulating con-stitutive relations for the momentum supply terms, which in fact may include concon-stitutive variables from the interacting constituents, the imposed restrictions on the momentum supply terms can be utilized to formulate an extended equation of motion for the con-stituents. Since no obvious way to constitute the interaction terms exists, eliminating these supply terms between the constituents simplifies the model considerably. The main advantage of such a formulation given here is that no new constitutive relations need to be formulated and that the existing calculation scheme developed by [3], can easily be supplemented and used.
The motion of the fluid within the foundation affects the strength and the deformation characteristics of the soil material. One of the first presented relations for characterizing the motion of the fluid is Darcy’s law, which is an empirical relation based on experi-mental study of steady flow in a vertical column of homogeneous sand. The extension to three dimensions and to unsteady flow were suggested as heuristic generalizations and verified from experiments. However, attempts have been made to derive Darcy’s law, or more generally, the equations of motions in porous media, from basic principles of hydro-dynamics using different conceptual model approaches. Among these conceptual models are the capillary tube models, the capillary fissure models, the hydraulic radius models, resistance to flow models, etc., c.f. [4]. A few attempts have been made to derive the equation of motion from Navier-Stokes equation. The Navier-Stokes equation has been solved making approximations or by averaging, e.g. [5], [6]. The different derivations of Darcy’s law resulted in different expressions for the permeability constant, depending on the geometry of the conceptual model used in the derivation. It has been also shown by experiments that the Darcy’s law holds strictly only for laminar flows.
equa-tion of the soil is derived including the acequa-tion of body forces, the stresses in the fluid and the effects of the fluid flowing. The magnitudes and directions of these terms are discussed and approximate expressions are derived. Although the mixture model con-sists of two constituents it is possible to analyze a system with three constituents. The constitutive model used for the soil constituent includes voids and proper scaling of the volume fraction enables one simulation of tests on water unsaturated soil. The structure of the model benefits from computer algorithms for integrating the constitutive relation for granular materials developed in [3], and after slight modifications it is used in finite element simulations of water saturated sand.
2 BINARY MIXTURE MODEL FOR SOIL PHYSICS
2.1 Kinematics and balance equations
In this section the kinematics and kinetics of a mixture of N constituents will be discussed briefly. The notations and definitions necessary for the implementation of the model derived from a mixture theory will be presented. If Bt denotes the region in
space occupied by the mixture at time t, each spatial position x at time t of the current placement of the mixture is simultaneously occupied by particles Xα from all of the
constituents α, α = 1, . . . , N . For each constituent the motion is given by x = χα(Xα, t).
The vector x is the spatial position occupied by the particle Xα at time t and Xα is the
position in the reference configuration. The functionχα is the deformation function of the α-th constituent and it is assumed to be invertible, i.e. Xα = χ−1α (xα, t). The velocity
and the acceleration for each of the constituents of the mixture at time t are defined by ´
xα = (∂χα(Xα, t)(/)∂t) , ˝xα = (∂2χα(Xα, t))/(∂2t)(1)
The deformation gradient for constituent α at time t is a linear transformation given by,
Fα = Gradα χα(Xα, t) = (∂χα(Xα, t))/(∂Xα) (2)
where the differential operator “Grad” denotes partial differentiation with respect to the reference position Xα. The inverse of the deformation gradient is given by
F−1α = gradχ−1α = (∂χ−1α (x, t))/(∂x) (3)
where the differential operator “grad” refers to partial differentiation with respect to the spatial point x.
In general, the individual constituents undergo different motions and, thus, material time derivatives, related to the motion of each constituent α, need to be introduced. If Γ(x, t) is an arbitrary differentiable function, the material time derivative, following the motion of the α-constituent, is defined by,
where the grad-operation on Γ denotes the partial derivative of Γ with respect to the local position x. The material velocity gradient of the α-constituent is then given by,
(Fα)α = Grad ´xα = (∂ ´xα)/(∂Xα) (5)
whereas the spatial velocity gradient is defined by Lα = grad ´xα. When formulating
constitutive equations the velocity gradients need to be decomposed into a symmetrical part, the rate of deformation tensor Dα, and a skew-symmetric part, the spin tensor
Wα, so that Lα¯Dα + Wα, with Dα = symLα = 12(Lα + LTα) and Wα = skewLα = 1
2(Lα− L T
α). For a mixture, each constituent is assigned a density and ρα will denote the
density of the α-th constituent. Then ρα is given at each position x ∈ Bt at any time t,
i.e. ρα = ρα(x, t). The density of the mixture, ρ, at x and at time t is then defined by
ρ = ρ(x, t) = Nα=1ρα(x, t).
The density ρα represents the mass dmα of the α-th constituent per unite volume of
the mixture, and ρα is referred to as the bulk, or partial, density. The true or material
density of the constituent, which is denoted by γα, is the mass per unit of the true volume
dVα of the α-th constituent, ρα = dmα/dV and γα = dmα/dVα. The volume dV of the
mixture is the sum of the partial volumes, i.e. dV =
N
α=1
dVα. The volume fraction φα
of the α-th constituent is defined as the local ratio of the volume element dVα of a given
constituent α with respect to the volume element dV of the mixture, φα(x, t) = dVα/dV
Obviously the volume fraction φα is related to the two density functions ρα and γα
through the relation φα(x, t) = ρα(x, t)/γα(x, t). If there are no voids in the mixture, the
saturation condition follows as Nα=1φα(x, t) = 1 The volume fraction and its spatial
and temporal changes may affect the response functions in the constitutive relations for an immiscible mixture. The concept incompressible for a mixture is not a straightforward notion. If the α-th constituent is incompressible then γα is a constant. The bulk density
ρα, however, need not be constant even if the α-th constituent is incompressible.
Finally, for the mixture the barycentric velocity, or simply the velocity of the mixture at (x, t), is introduced as the weighted average of the constituent velocities defined by
˙x = ˙x(x, t) = 1 ρ N α=1 ραx´α(x, t) (6)
The acceleration of the mixture is defined by ¨x = ¨x(x, t) = ∂t˙x + L ˙x. and the
diffusion velocity of the α-constituent is defined as wα(x, t) = ´xα(x, t) − ˙x(x, t). In
Conservation of mass ∂ρS ∂t + divρS´xS = 0 ∂ρF ∂t + divρFx´F = 0 (7) Balance of momentum ρF x˝F = div TF + ρF bF + ˆpF ρS x˝S = div TS + ρS bS + ˆpS (8) Balance of energy ρFε´F = trTFLF − divhF + ρFrF + ˆεF ρSε´S = trTSLS − divhS+ ρSrS + ˆεS (9)
TS and TF are the Cauchy stress tensors for the constituents and here it is assumed
that the partial stress tensors are symmetric, that is TST=TS and TFT=TF. bS and
bF are the specific body forces, and ˆpS and ˆpF are the momentum supply terms due to
interactions between the constituents. In the energy balance equation εS and εF are the
specific internal energies, hS and hF the partial heat flux vectors due to conduction, rS
and rS are the specific influxes of heat due to radiation, and ˆεS and ˆεS the energy supply
terms due to interactions between the constituents. The corresponding set of balance equations for the mixture reads:
∂ρ
∂t + divρ ˙x = 0 (10)
ρ¨x = divT + ρb (11)
ρ ˙ε = trTL − divh + ρr (12)
where T is the Cauchy stress tensor, b the specific body force, ε the specific internal energy, h the heat conduction vector and r the influx of heat from radiation of the mixture.
The balance of energy equation will not be treated in this work and in the following we concentrated only on the balance of momentum. By adding (7)-(8) for the two constituents and compare with equations (10)-(11) for the mixture, the following restrictions and definitions, relating the total and partial quantities are obtained,
ˆ
pS + ˆpF = 0 (13)
b = 1
ρ (ρS bS + ρF bF) (14)
T = TS + TF − ρSwS ⊗ wS − ρFwF ⊗ wF (15)
where wSand wF are the diffusion velocities of the soil and fluid constituents, respectively
2.2 Equation of motion for the soil
The acceleration fields for the fluid, ˝xF, and the soil, ˝xS, can be expressed by use of the
conservation of mass (7). It will be assumed that both the soil and the fluid motions are in quasi-steady state. For the soil it is further assumed that the velocity ´xS is negligible,
i.e. ´xS 0 . The approximated acceleration fields for both constituents then become,
ρF˝xF = ∂(ρF∂tx´F) + div(ρF x´F ⊗ ´xF) div(ρF x´F ⊗ ´xF) ρSx˝S = ∂(ρ∂tSx´S) + div(ρS x´S⊗ ´xS) 0 (16)
Then the balance of momentum (8) can be approximated by div(ρF x´F ⊗ ´xF) = div TF + ρF bF + ˆpF 0 = div TS + ρS bS + ˆpS (17)
Combining the two equations in (17) and using the restriction relations (13) and (14), enables elimination the momentum supply terms, giving the following relation,
0 = div (TS + TF − ρF x´F ⊗ ´xF) + ρ b (18)
To constitute the stress tensor for the fluid the classical Navier-Stokes relation is employed, whereas to constitute the stress tensor for the soil the elasto-plastic theory developed in [2] and briefly described in Section 2.3 is applied.
Compared to other mixture models, e.g. [7], the advantage of the last equation is that the supply term does not need to be constituted and that in the numerical simulation the existing scheme for integration of the stresses for dry sand can be easily supplemented and employed. Having computed the stress tensor TS for the soil, it can be corrected
by the last three terms in (18) to obtain a new stress tensor for the soil including the influence of the fluid.
A similar procedure as for the balance of momentum equation can be applied to the balance of energy equation. This will result in an energy equation, where the supply terms are eliminated.
2.3 Constitutive relations
In order to have a closed system of equations the postulated conservation and balance equations have to be supplemented with a set of constitutive equations. As the set of dependent constitutive variables, the partial stress tensors, TS, TF, and the momentum
the fluid is constituted as an incompressible linearly viscous fluid and the stress tensor
TF takes form:
TF = −pF 1 + 2µDF (19)
DF is the rate of deformation tensor of the fluid and the viscosity coefficient µ is a material
parameter.
In the present work, the non-associated small deformation plasticity model developed in [2]. The model is formulated for small deformations and without including rate effects. The material model is calibrated and validated to experimental data for dry sand in [8]. The small deformation model showed good capability to represent the mechanical behaviour of sand, suggesting that reliable results could be obtained also in the present study. The limited number of material parameters, six only, together with the numerical algorithm developed in [3], makes the model very useful for solving geotechnical problems of engineering interest using the finite element method. In the present study the rate effects are not included, but if the material model for the soil is reformulated or if another material model for soil able to handle rate effects is incorporated in the mixture model, then the proposed method will be able the represent loading rate effects. The non-associated small deformation plasticity model have been generalized to a large deformation formulation in e.g. [10] and the model has been formulated using the Jaumann rate of the Cauchy stress tensor for the soil. A brief description of the constitutive model for small deformation is given below. The Cauchy stress tensor TS of the soil is formulated using
the stress invariants. The model operates with the traditional split of the stress tensor
TS and the strain tensor ES into hydrostatic and deviatoric parts
pS =−1 3I1 =− 1 3trTS, SS = TS+ pS1, EvS =−trES, E˜S = ES+ 1 3EvS1 (20) where pS is the mean pressure, I1 is the first stress invariant, SS the deviatoric stress
tensor, EvS volumetric strain and ˜ES the deviatoric strain tensor. Restricting to small
deformations at each point, the displacement field u(x, t) and its gradient remain small as compared to the macroscopic reference length of the model. Using this geometrically linear approach the strain tensor ES is additively split into an elastic, EeS, and a plastic,
EpS, part,
ES = EeS + EpS (21)
In order to describe the plastic material properties the yield function f (TS) is needed to
bound the elastic domain. In terms of the mean pressure pS and the third stress invariant
I3 = det TS, the yield function in the model is expressed as,
f (TS) = I3+ p3S(pS/pf S)m (22)
The yield function grows in self-similar way. The size of the yield function is controlled by the parameter pf S, which is the only hardening parameter in the model, and the exponent
in the deviatoric plane, whereas the volumetric part is non-associated, leading to a similar format as for the yield surface,
g(TS) = I3+ p3S(1− γg2(pS)) , γg(pS) = 1− (pS/pgS)n (23)
where γg is a shape function and the exponent n is a material constant.
The elastic constitutive tensor Ce has the same format as in the Cam-Clay model, cf. [9], and is written in tensor form as,
Ce = (K − (2/3)G)1 ⊗ 1 + 2GI (24)
where K = pS/κ and 1 and I are the second and fourth order unit tensors, respectively.
The shear modulus G and the non-dimensional flexibility parameter κ are material con-stants.
The direction of the plastic strain increment dEpS is the gradient of the plastic potential
function and its magnitude is given by the plastic multiplier dΛ
dEpS = dΛ ∂g/∂TS (25)
The hardening of the yield surface depends on both plastic volumetric and deviatoric strain increments in order to incorporate the possibility to model dilatancy. Thus, the hardening rule in the model is a weighted sum of the volumetric and the deviatoric parts of the plastic work,
dpf S = 1 λ− κ pSdEpvS + w SS· d˜EpS (26) where w and λ are material constants. The hardening of the yield function in the present model is controlled by the size parameter pf S, and thus, the hardening parameter H =
∂f /∂Λ becomes
H = − (∂f/∂pf S) (∂pf S/∂Λ) (27)
The elasto-plastic stiffness matrix Cep is needed for use in the global, non-linear equation system and reads
Cep = Ce − (H + ∂f/∂TS· Ce ∂g/∂TS)−1(Ce ∂g/∂TS⊗ ∂f/∂TS Ce) (28)
The Cauchy stress tensor for the soil TS is there after obtained from
TS = CepES (29)
3 EVALUATION OF THE BODY FORCES, THE RATE OF
DEFORMA-TION TENSOR AND THE VELOCITY OF THE FLUID
After insertion of the constitutive relations for the fluid and decomposition into devi-atoric and hydrostatic parts, relation (18) reads,
where p = pS+ pF is the total hydrostatic pressure. In order to compute the total stress,
a suitably way of estimating the rate of deformation tensor DF, the velocity field ´xF and
the body force b is needed.
By assuming that there exists a potential Φ to the specific body force, i.e. b = grad Φ, the last term in (30) is expressed as,
ρ b = ρ grad Φ = grad(ρΦ) − Φ grad(ρ) (31)
If the changes in the density ρ over the body are assumed to be small, i.e. gradρ 0, we obtain,
ρ b grad(ρΦ) = div(ρΦ 1) (32)
The approximated expressions inserted into (30) yields,
0 = div (SS − (p − ρΦ) 1 + 2µDF − ρF x´F ⊗ ´xF) (33)
The body force acting on the fluid is the gravity force and, thus, Φ =−gz, where g is the gravity acceleration and z is the coordinate in the direction of g.
In order to estimate the last term in (33), both the magnitude and the direction of the velocity of the fluid must be determined. A possible way of estimating the magnitude of the velocity of the fluid is through an approximated equation of balance of momentum of the fluid constituent. To obtain an estimate of the magnitude of the velocity field of the fluid flow, the equation of balance of momentum for the fluid is approximated by
¯
ρF x˝F = −grad pF (34)
where ¯ρF, the true density of the fluid, is introduced and defined by ¯ρF = ρF/φF. The
so-called true balance of momentum (34) is appropriate for immicsible mixtures, see [11] and [12], and states that the motion merely is driven by the gradient of the partial pressure pF, neglecting the viscoeffects. This is an approximate equation of balance of momentum,
in which the deviatoric part of the stress and the interaction term are ignored.
Assuming quasi-steady motion and the absence of vorticity, ωF, the acceleration of the
fluid constituent will be ˝ xF = ∂ ´xF ∂t + LFx´F = ∂ ´xF ∂t + 1 2 grad(´xF) 2+ 2(ω F x ´xF) 1 2 grad(´xF) 2 (35)
The approximations in the last equation are made because the process investigated here is monotonic or stationary and without turbulence. However, the intermediate expression in (35) is general and can be used if turbulent and non-stationary flow is considered.
If UF denotes the speed of the fluid constituent, combining (34) and (35) yields
In this equation it has been assumed that the changes in ρF are small. This is also
confirmed by the simulations presented in the last section.
Integrating (36) together with the initial condition pF = pF 0at t = 0 gives an estimate
of the speed UF of the fluid as
UF2 = 2
¯
ρF (pF 0− pF) (37)
The partial pressures of the fluid constituent pF and pF 0 can be related to the total
pressure and volume fractions through the relations pF = φFp and pF 0 = φF 0p, where
φF 0 is the volume fraction of the fluid at the initial state. Insertion into (37) yields an
expression for UF as function of the volume fraction of the fluid φF:
UF2 = 2 p ¯
ρF (φF 0− φF) (38)
The saturation condition φF + φS = 1 allows for writing UF as a function of the volume
fraction of the soil as
UF2 = 2 p ¯
ρF (φS− φS0) (39)
This result can also be derived from the Bernoulli equation by considering a steady flow of an incompressible fluid with no gravity, cf. [13].
Employing (39), the magnitude of the term ρF x´F ⊗ ´xF in (33) can be estimated:
|ρF x´F ⊗ ´xF| ρFUF2 2p(φS− φS0) (40)
In addition, the direction of the velocity needs to be known. There is no obvious way of how to choose the direction of the velocity field of the fluid. The aim here is to simplify the mixture model in order to facilitate use in simulations of engineering interest.
For that reason the velocity of the fluid is chosen to have its change of direction related to the principal strain increment dF of the fluid at each load step, which consecutively
is approximated by the principal strain increment dS of the soil. An unit vector n is
computed as,
n = d S1e1+ d S2e2+ d S3e3
(d S 12+ d
S 32+ d S 32)1/2
(41) where e1, e2 and e3 are the principal strain directions. Then the velocity of the fluid is
estimated as,
´
xF − UF n (42)
Thus it is assumed that the streamlines are parallel to the principal strain increment but in opposite direction.
The remaining term in (33) includes the rate of deformation tensor DF, defined as the
symmetric part of the spatial velocity gradient LF, i.e.:
DF = symLF = 1
2(LF + L
T
There are two possible ways of evaluating LF. One can proceed from the velocity
estimated in (42) and calculate LF. Using the deformation gradient FF, the spatial
velocity gradient for the fluid LF is, from (??), defined by,
LF = ˙FF F−1F (44)
In the numerical simulations the deformation gradient is calculated at the beginning and the end of each load increment so that, for a time increment ∆t = t2 − t1, the time
derivative ˙FF is approximated by,
˙
FF = (FF(t2) − FF(t1))/(∆t) (45)
4 NUMERICAL ALGORITHM
The present binary mixture model formulation facilitates the use of the finite element technique, which is the most practical for strength and deformation analysis of soil bodies. The formulation of the mixture model through deriving the equation of motion is that an existing numerical algorithm, developed in [3], for integration of the stress tensor for dry sand can be used by making only minor modifications. The implicit integration procedure implemented as a user defined mechanical material behaviour in the finite element code Abaqus is utilized. The stress tensor and the hardening variable are determined through solving a non-linear system of equations iteratively. Once the stress tensor for the dry sand is integrated, it is supplied by the terms in (33) and the tangent stiffness matrix is updated. A brief description of the numerical algorithm used to integrate the stress for the sand is given below.
The constitutive calculations are performed using the implicit integration algorithm for-mulated for models for granular materials in [3]. In the implicit backward-Euler method the solution makes use of the gradient to the potential surface at the final stress state, which is unknown. The final stress and hardening parameters are determined solving the non-linear equations iteratively, so that the stress increment fulfills the consistency condition. The integration scheme determines the stress changes ∆TS and the hardening parameters corresponding to a total change of displacement ∆ES within the load incre-ment. In the implicit backward-Euler method, the solution makes use of the gradient to the potential surface at the final stress state (TSn, pfS n), which is unknown. The
cur-rent estimated state, (TS, pfS), is determined in each iteration step relative to the former
equilibrium state (TS n−1, pfSn−1). The computed stress and hardening parameters must fulfill the consistency relation,
f (TS n, dpfSn) = 0 (46)
The total strain increment, ∆ES = ES n− ES n−1, is the sum of the elastic and plastic
parts,
The elastic part is written as ∆ESe =ESe(TS n)− ESe(TS n−1), whereas for the plastic
part use of a non-associated flow rule implies ∆ESp = ∆χ [∂g(TS n)/∂TS]. Insertion into
Eq. (47) yields
ES n− ES n−1 =ESe(TS n)− ESe(TSn−1) + ∆χ (∂g(TSn))/(∂TS)(48)
The stresses TS n and the plastic multiplier χ are unknowns, whereas the prescribed total strain ES n and terms related to the previous equilibrium state (n− 1), are known. Since isotropic hardening is assumed pfS is the only hardening parameter in the model,
controlling the size of the yield function, and it must be updated continuously,
pfS n− pfS n−1 = H2(TS n) ∆χ (49)
The non-linear equation system comprising Eqs.(46), ( 48) and (49) can be solved using the Newton-Raphson iteration scheme. By a first order Taylor expansion about the current state (TS n, pfSn), these three equations are linearlized. The hardening parameter pfS
can be calculated explicitly at each intermediate state, which allows for a reduction of the number of equations in the iteration scheme and we obtain the following two linear equations that have to be solved,
Ce−1 + ∂ 2g ∂TS∂TS δTS + ∂g ∂TS∂TS δχ = δES ∂f ∂ES − ∆χ ∂H2 ∂TS δTS − H1H2δχ = −f(TS, pfS) (50)
where ∂ESe/∂TS = Ce−1 have been used, with Ce denoting the elastic constitutive
matrix as defined in Eq.(24). After solving this linear system of equations in terms of (δTS, δχ), the hardening parameter δpfS is updated explicitly at each iteration. Finally,
the increments (∆TS, ∆χ) are updated by subincrements until the convergence is reached,
∆TS = ∆TS + δTS , ∆χ = ∆χ + δχ (51)
5 NUMERICAL SIMULATIONS AND DISCUSSION
no voids, i.e. a soil system containing soil and water only. In the numerical simulations a monotonic loading process is enforced. The simulations showed that the influence of the correction term 2µDF is insignificant, and this term is omitted in the numerical examples
presented in next section. The influence of the gravity force is of varying significance, depending on the specific test simulated. Since the body force depends of the depth, near the surface it may be insignificant, whereas for an analysis of soil placed deep beneath the surface, it has a considerable influence. The following examples concentrate only on the influence from pore water pressure and flow of the water. The capability of the mixture model is illustrated at two levels. The first example concerns a triaxial test on one element and the second example is a multi element test of a flexible footing resting on the surface of saturated sand.
5.1 One element test
A triaxial test on one linear brick element has been conducted. The boundary con-ditions for the soil and the applied load are shown in Fig. 1. The test starts at initial hydrostatic pressure p0 = 800 kPa. This initial condition was imposed on the three loaded
surfaces. Additionally, an uniformly distributed load of 2500 kPa was applied on the top surface in the vertical direction, as shown in Fig. 1. The following material constants
Initial surface load
Initial surface load Additional surface load
Initial surface load
Figure 1: Boundary and load conditions for triaxial test on one element.
are relevant for loose sand with void ratio e = 0.85: G = 47300 kPa, λ = 0.00759, κ = 0.00609, n = 1.0793, m = 0.4498 and w = 0.3066, cf. [8]. In this study these parameters were used for simulations of water unsaturated sand. When analyzing water saturated sand, κ and λ were scaled by the specific volume v = 1.85 to λ = 0.00655 and κ = 0.00327. The initial volume fraction of the sand in the water saturated analysis becomes φS = 0.54, given the void ratio e = 0.85. The boundary conditions for the water
At the beginning At the end of the test of the test
Water 46% 47.6%
Voids 25% 24%
Sand 29% 28.4%
Saturation 65% 66.3%
Table 1: Initial and final saturation for the water unsaturated example.
respect to the strain of the sand according to Eq.( 42), the water is allowed to flow into the volume through the front and right surfaces in Fig. 1, whereas the water is allowed to flow out through the top surface. The volume fraction of the sand changes during the test as the sand compacts in the beginning of the test, whereafter occurring at an axial strain of 0.62, starts to dilate, see Fig. 4.2b. Therefore it is assumed that the water flows into the element equally through the front and the right boundary, when the sand is compacted, whereas it flows out through the top boundary after the sand has started to dilate.
For the example with unsaturated soil, the degree of saturation is included in the in the stiffness parameters κ and λ. κ and λ are scaled by the specific volume v = 1.85, which sub-sequentially determines the degree of saturation since the volume fraction of the soil is 0.54 in the example, i.e. the initial volume fraction of the water is 0.46. In the presented example for the water unsaturated soil, we have 71 % voids whereof 46% are filled with water, yielding a degree of saturation of 65 %. The degree of saturation varies slightly during the test example, and at the end of the simulation it is 66.3%. This is due to compaction of the sand during the test, i.e. the volume fraction of the sand has decreased slightly. The void ratio of 0.85 is kept during the test, resulting in 24% voids at the end of the test, see Table 5.1. The total volume fraction, i.e. the sum of the volume fraction for soil, water and voids, is kept constant, φS+φF+φV = 1, which leads
to a slight increase in volume fraction of water, i.e. in saturation. The initial saturation is 65% and it increases slightly during the test to 66.3%. At the end of the test, some of the voids have been filled with water. For the water unsaturated example the boundary condition for the water are the same as in the simulation with water saturated sand.
The water unsaturated soil can be considered as consisting of three constituents. Al-though the model was formulated as a binary mixture model, it is still possible to analyze water unsaturated soil because the ratio between the voids and the soil was kept constant throughout the simulation.
the water unsaturated sand reaches a lower ultimate state as compared to the dry sand. This is due to the increase in the specific volume v by adding water to the system. The sand becomes looser and, as expected, reaches a lower ultimate state. In comparison of the water saturated sand to dry sand it is seen that the curve for the water saturated sand is steeper before reaching the same ultimate state. This results is expected, since a mixture containing sand and water is stiffer than dry sand. Comparing the volumetric curves between water unsaturated and dry sand, it is seen in Fig. 2b) that they coincide at the beginning of the test, but subsequently the volumetric strain becomes higher in the case of water unsaturated sand. The water between the sand particles results in a higher degree of compaction of the sand. This can be explained partly by the increase in specific volume v by adding water to the system. The sand becomes looser, which makes it easier to compact. Also water decreases the apparent friction between particles, which allows them to compact more in the presence of voids of air. The volumetric curve of water saturated sand is lower than that for the dry sand, i.e. the water between the sand particles results in a lower degree of compaction of the sand. Neglecting the body forces and viscoeffects, the equation of motion (33) becommes,
0 = div (TS − pF1 − ρF x´F ⊗ ´xF) (52)
where the total hydrostatic pressure p has been separated into partial hydrostatic pressures pS and pF. The difference between the total stress and hydrostatic pressure of the fluid pF
is the so-called effective stress, usually used in soil mechanics, cf. [14], and the response of the soil depends on the effective stress. The interpretation of (52) is that the total stress contains an additional term, −ρF x´F ⊗ ´xF, accounting for the fluid flow.
Figs. 3 and 4 show results for water saturated sand. In Fig. 3, the effects of the pore water pressure and the flow of the water are compared. The stress-strain curves in Fig. 3a) are, in practice, identical with an unsignificant influence on the ultimate state. There is, however, observable differences between the volumetric curves. The presence of flowing water results in a higher degree of compaction, which is expected since the movement of water makes it easier for the soil to compact. In Fig. 4a) the hydrostatic pressure of the
0 0.025 0.05 0.075 0.1 0 0.2 0.4 0.6 0.8 Vertical strain ε 3 Stress ratio ( σ 3 − σ 1 )/ σ 3 (MPa)
Water saturated sand Dry sand
Water unsaturated sand
0 0.025 0.05 0.075 0.1 −5 0 5 10 15x 10 −3 Vertical strain ε 3 Volumetric strain ε v = ε 1 + ε 2 + ε 3
Water saturated sand Dry sand
Water unsaturated sand
0 0.025 0.05 0.075 0.1 0 0.2 0.4 0.6 0.8 Vertical strain ε 3 Stress ratio ( σ 3 − σ 1 )/ σ 3 (MPa)
including pore water pressure
including pore water pressure and water movement
0 0.025 0.05 0.075 0.1 −2 0 2 4 6 8x 10 Vertical strain ε 3 Volumetric strain ε v = ε 1 + ε 2 + ε 3
including pore water pressure including pore water pressure and water movement
Figure 3: a) Stress ratio and b) volumetric water saturated sand, versus vertical strain .
0 0.025 0.05 0.075 0.1 0.5 0.75 1 1.25 1.5 Vertical strain ε 3 Pressure (MPa)
Pore water pressure Total pressure 0 0.025 0.05 0.075 0.1 0.53 0.535 0.54 Vertical strain ε 3
Volum fraction of sand
φ s In water saturated sand
In dry sand
Figure 4: a) Total pressure and pore water pressure, b) comparison of volume fraction of sand, versus vertical strain.
water saturated sand and the pore water pressure are visualized, whereas Fig. 4b) shows the changes of volume fraction of sand for the tests on dry and water saturated sand.
5.2 Multi element test
The simulation of multi element test is a tradition geotechnical problem of a flexible footing resting on a surface of water saturated sand. The symmetry allows half of the geometry only to be modelled. The right boundary in Figs. 5-10 is a reflected boundary, whereas the left and the bottom boundaries are smooth. The water is allowed to flow through the top boundary of the footing only. It is in the simulation assumed that the total volume fraction, i.e. the sum of the volume fraction for soil and water, is constant, i.e. φS+φF = 1, and that the direction of the water flow is determined by Eq.(42). A
S22 VALUE -2.18E+00 -2.07E+00 -1.96E+00 -1.86E+00 -1.75E+00 -1.64E+00 -1.53E+00 -1.43E+00 -1.32E+00 -1.21E+00 -1.10E+00 -9.97E-01 -8.90E-01 -7.82E-01 S22 VALUE -2.51E+00 -2.38E+00 -2.24E+00 -2.11E+00 -1.98E+00 -1.84E+00 -1.71E+00 -1.58E+00 -1.44E+00 -1.31E+00 -1.18E+00 -1.04E+00 -9.11E-01 -7.78E-01
Figure 5: Vertical stress distribution of footing test on dry sand (top) and water saturated sand (bottom).
S11 VALUE -1.90E+00 -1.80E+00 -1.70E+00 -1.59E+00 -1.49E+00 -1.39E+00 -1.28E+00 -1.18E+00 -1.08E+00 -9.73E-01 -8.70E-01 -7.67E-01 -6.63E-01 -5.60E-01 S11 VALUE -2.50E+00 -2.35E+00 -2.19E+00 -2.04E+00 -1.89E+00 -1.74E+00 -1.58E+00 -1.43E+00 -1.28E+00 -1.13E+00 -9.75E-01 -8.22E-01 -6.70E-01 -5.17E-01
Figure 6: Horizontal stress distribution of footing test on dry sand (top) and water saturated sand (bottom).
the normal stress in the three coordinate directions are shown. Besides the differences in pattern, all maximum values for the water saturated sand are higher than for dry sand. The largest difference, about 30%, in maximum value is observed in the horizontal stress value shown in Fig. 7. The displacement in the vertical direction is presented in Fig. 8. The two simulations follow the same pattern, but as reads the displacement beneath the footing, it is about 17% lower for the case of water saturated sand as compared to dry sand.
dis-S33 VALUE -1.61E+00 -1.55E+00 -1.49E+00 -1.42E+00 -1.36E+00 -1.29E+00 -1.23E+00 -1.16E+00 -1.10E+00 -1.04E+00 -9.72E-01 -9.08E-01 -8.43E-01 -7.79E-01 S33 VALUE -1.98E+00 -1.88E+00 -1.79E+00 -1.70E+00 -1.61E+00 -1.51E+00 -1.42E+00 -1.33E+00 -1.24E+00 -1.14E+00 -1.05E+00 -9.58E-01 -8.65E-01 -7.72E-01
Figure 7: Normal to the vertical stress distribution of footing test on dry sand (top) and water saturated sand (bottom). U2 VALUE -6.07E-02 -5.60E-02 -5.13E-02 -4.66E-02 -4.19E-02 -3.72E-02 -3.25E-02 -2.78E-02 -2.31E-02 -1.84E-02 -1.37E-02 -8.95E-03 -4.24E-03 +4.63E-04 U2 VALUE -5.03E-02 -4.64E-02 -4.24E-02 -3.85E-02 -3.45E-02 -3.06E-02 -2.66E-02 -2.27E-02 -1.87E-02 -1.47E-02 -1.08E-02 -6.84E-03 -2.89E-03 +1.07E-03
Figure 8: Vertical displacement of footing test on dry sand (top) and water saturated sand (bottom).
tribution of mean stress of the footing problem on dry sand at vertical load 1950 kPa is presented. As seen from Fig. 9 and Fig. 10, the overall pattern is the same, but for the dry sand the pattern appears at about 22% higher load.
6 CONCLUSIONS
PRESS VALUE +7.44E-01 +8.31E-01 +9.18E-01 +1.01E+00 +1.09E+00 +1.18E+00 +1.27E+00 +1.35E+00 +1.44E+00 +1.53E+00 +1.62E+00 +1.70E+00 +1.79E+00 +1.88E+00 PRESS VALUE +7.62E-01 +8.57E-01 +9.52E-01 +1.05E+00 +1.14E+00 +1.24E+00 +1.33E+00 +1.43E+00 +1.52E+00 +1.62E+00 +1.71E+00 +1.81E+00 +1.90E+00 +2.00E+00
Figure 9: Mean stress distribution of footing test on dry sand (top) and water saturated sand (bottom).
PRESS VALUE +7.10E-01 +8.35E-01 +9.60E-01 +1.08E+00 +1.21E+00 +1.33E+00 +1.46E+00 +1.58E+00 +1.71E+00 +1.83E+00 +1.96E+00 +2.08E+00 +2.21E+00 +2.33E+00
Figure 10: Mean stress distribution of footing test on dry sand at load 1950 kPa.
also been utilized to simulate a multi element boundary value problem. The traditional geotechnical problem of footing resting on a soil surface have been simulated and the results has been compared to simulations on dry sand.
REFERENCES
[1] A.Ahadi and E.Lundgren,A binary mixture model for soil physics - a set of novel governing equations, Report, Mechanics, Lund University, Sweden,(2003).
[2] S. Krenk. Characteristic State Plasticity for Granular Materials, Part1: Basic The-ory. Int J Solids and Struct, 37, 6343-6360, (2000).
[3] A. Ahadi and S. Krenk. Implicit Integration of Plasticity Models for Granular Ma-terials, Comp. Methods. Appl. Mech. Engng., 192, 3471–3488, (2003).
[4] J.Bear.Dynamics of Fluids in Porous Media, Dover Publications,New York,(1972). [5] S. Whitaker. The equation of motion in porous media,Chem. Eng. Sci., 21, 291–
300,(1972).
[6] S. Irmay. On the theoretical derivation of the Darcy and Forchheimer formulas, . Amer. Geophys. Union No. 4, 39, 702–707, (1958).
[7] W. Ehlers and T. Graf. Deformation and localization analysis of unsaturated soil, Proc. 5th World Congress on Comput. Mechanics (WCCM V), 2002, Vienna, Aus-tria, Editors: Mang, Rammerstorfer, Eberhardsteiner,(2002).
[8] A. Ahadi and S. Krenk. Characteristic State Plasticity for Granular Materials, Part2: Model Calibration and Results. Int J Solids and Struct, 37, 6361-6380, (2000). [9] A.N.Schofield and C.P. Wroth. Critical state soil mechanics, McGraw-Hill, New
York, (1968).
[10] O. Kristensson and A. Ahadi. Numerical study of localization in soil sys-tems,Computers and Geotechnics, 32, 600-612, (2005).
[11] A.Bedford and D.S. Drumheller. Theories of immiscible and structured mixtures. Int. J. Engng. Sci., 21, 863-960,(1983).
[12] D.S. Drumheller. On theories for reacting immiscible mixtures.Int J of Eng Scin,
38,347-382,(2000).
[13] G.E.Mase and G.Mase. Continuum Mechanics for Engineers, CRC Press,2nd,(1999). [14] D.M. Wood. Soil Behaviour and Critical State Soil Mechanics, Cambridge University
