Next Article in Journal
Geo-Heritage Specific Visibility as an Important Parameter in Geo-Tourism Resource Evaluation
Previous Article in Journal
Performance of Remotely Sensed Soil Moisture for Temporal and Spatial Analysis of Rainfall over São Francisco River Basin, Brazil
Previous Article in Special Issue
Nanoindentation Studies of Plasticity and Dislocation Creep in Halite
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Finite Element Simulations of an Elasto-Viscoplastic Model for Clay

by
Mohammad N. Islam
1,*,
Carthigesu T. Gnanendran
1 and
Mehrdad Massoudi
2
1
School of Engineering and Information Technology, University of New South Wales, Canberra, ACT 2612, Australia
2
US Department of Energy, National Energy Technology Laboratory (NETL), P.O. BOX 10940, 626 Cochrans Mill Road, Pittsburgh, PA 15236, USA
*
Author to whom correspondence should be addressed.
Geosciences 2019, 9(3), 145; https://doi.org/10.3390/geosciences9030145
Submission received: 19 February 2019 / Revised: 14 March 2019 / Accepted: 19 March 2019 / Published: 26 March 2019
(This article belongs to the Special Issue Micromechanics of Reservoir and Cap Rocks)

Abstract

:
In this paper, we develop an elasto-viscoplastic (EVP) model for clay using the non-associated flow rule. This is accomplished by using a modified form of the Perzyna’s overstressed EVP theory, the critical state soil mechanics, and the multi-surface theory. The new model includes six parameters, five of which are identical to those in the critical state soil mechanics model. The other parameter is the generalized nonlinear secondary compression index. The EVP model was implemented in a nonlinear coupled consolidated code using a finite-element numerical algorithm (AFENA). We then tested the model for different clays, such as the Osaka clay, the San Francisco Bay Mud clay, the Kaolin clay, and the Hong Kong Marine Deposit clay. The numerical results show good agreement with the experimental data.

1. Introduction

Saito and Uezawa Saito and Uezawa [1] showed that the Takabayama landslide in Japan was mainly due to the time dependent deformation of the clay, resulting in the failure of the slope. In 1963, the failure of the Vayont reservoir in Italy, cost more than 3000 lives [2]; this was again due primarily to the creep of clay. These and many other examples in geotechnical engineering, reservoir geomechanics, wellbore plugging in petroleum and mining industry, etc., have made studying and modeling of clay a very important topic.
Clay is considered to be a multi-component system composed of solid particles, and fluid components (gas and/or liquid) [3]. For many contaminated soils, there might be additional components, such as bubbles, oil spill, etc. In a partially saturated clay, the void space is filled with gas and/or liquid, and, in a fully saturated clay, only a liquid phase is present. Removal of any fluid from this system depends on the magnitude of the external load and the drainage path, the imposed duration, and the properties of the clay. When an external load is applied, pressure is developed in the fluid, which may dissipate instantaneously or gradually [4]. On the other hand, for a solid medium with non-crushable fine materials, when the fluid pressure dissipation is negligible, the settlement of clay does not end. In this case, creep may continue for a long time under a constant pressure [5]. This process constitutes a complex hydro-mechanical phenomenon [6,7].
Many constitutive models have been developed for clay. These range from simple elastic models to very complicated elasto-viscoplastic models. An overview can be found in the works of Owen and Hinton [8] and Desai and Siriwardane [9]. Chaboche [10] also presented a review of constitutive theories for plastic and viscoplastic type models (see also Tatsuoka et al. [11]). Liingaard et al. [12] suggested that the time-dependent viscous type models can be separated into: (i) empirical models (Singh and Mitchell [13]); (ii) rheological models (Feda [14]); and (iii) models considering “stress-strain-time” (Adachi and Okano [15]). The first two categories are important for the basic understanding of the behavior of clay. In clay-based research and its finite element implementation, the third category is most popular. For the remainder of this paper, the discussion is limited to the elasto-viscoplastic (EVP) type models. In the formulation of the EVP models, the strain rate consists of an elastic component and an inelastic component. To obtain an expression for the latter, it is essential to include the time-dependent viscous property in the model formulation. Sekiguchi [16] subdivided EVP models into: (a) the overstressed type [15]; (b) the non-stationary flow surface type (Sekiguchi [17]); and (c) others, such as the Bounding surface model [18] and Borja model [19]. Details of EVP models and their strength, as well as their limitations, can be found in the work by Liingaard et al. [12].
In this paper, we develop a new EVP model with a non-associated flow rule, considering Perzyna’s viscoplastic theory [20] and the MCC framework [21]. The time-dependent viscous aspect is incorporated using the Borja and Kavazanjian [19] approach. The new model requires six parameters. The current study is a generalization of the previous work by Islam and Gnanendran [22], who showed the strengths and the limitations of the existing methods when using the two-surface approach in the EVP models by incorporating the associated flow rule. Here, we use the non-associated flow rule and the three-surface approach. We also look at the Osaka clay, the San Francisco Bay Mud clay, the Kaolin clay, and the Hong Kong Marine Deposit clay in a triaxial loading situation. The wide range of experiments included the drained and the undrained tests, the strain rate test, the relaxation test, and the over consolidation effect tests. The model was implemented in a finite-element code using the numerical algorithm (AFENA) [23]. The pertinent details of the finite element implementation and the development of the non-associated flow rule are discussed in this paper.

2. Basic Equations

We consider clay to be a mixture of two constituents composed of solid particles and a liquid. We use the basic ideas and principles of continuum theories of mixtures, where electromagnetic and thermo-chemical effects are ignored (see Atkin and Craine [24] and Bowen [25]). The motion for each component is given through a mapping
x = χ i ( X i , t ) i = s , f
where χ i   designates the mapping, and i refers to s (solid component-clay) and f (fluid component). For clay, the kinematical quantities of interests are the displacement vector u , the velocity vector v s , the deformation gradient tensor F , and the linearized strain ε s :
u ( x , t ) = x ( x , t ) X s
v s = X s t
F = x s X s
ε s = 1 2 [ ( u x ) + ( u x ) T ]
where superscript T designates the transpose of a tensor. If V 0 is the volume of the particles before the deformation and V is the volume after the deformation, then a quantity of interest is the volumetric strain, also called the dilatation
ε v = t r   ε = ε k k = V V 0 V 0
Furthermore, in elastoplastic theories, for small deformations, we use ε ˙ instead of the symmetric part of the velocity gradient D s . That is, for small deformation, we assume ε ˙ = D s , where
ε ˙ i j = 1 2 ( u ˙ i x j + u ˙ j x i )
For the fluid component, the kinematical quantities of interest are the velocity vector v f , the gradient of the velocity L , and its symmetric part D , which are given by, respectively,
v f = X f t
L = v f x
D = 1 2 [ L + L T ]
We assume that the density of the solid and fluid components in their current configuration are ρ s and ρ f , respectively, while, in their pure (reference) state, i.e., before mixing, they are given as ρ R s and ρ R f . The density and the velocity of the mixture are then defined as
ρ ( x , t ) = ρ s ( x , t ) + ρ f ( x , t )
v ( x , t ) = 1 ρ [ ρ s v s + ρ f v f ]
Note that Equation (12) is one of the many possibilities for defining a mixture velocity. We can also define a relative velocity w such that
w ( x , t ) = v f ( x , t ) v s ( x , t )
The densities in the current and reference configuration are related through the kinematical field φ ( x , t ) , called the porosity, such that
ρ s = ( 1 φ ) ρ R s
ρ f = φ ρ R f
where 0 φ ( x , t ) φ m a x < 1 . Thus, when φ = 1 , we have a fluid with no pores and no particles. Note that this automatically assumes that the mixture is saturated. Otherwise, the porosities are constrained by φ s + φ f 1 . In geomechanics-related problems, we sometimes use the void ratio e , which is related to the porosity φ , through
e = φ 1 φ
Assuming no inter-conversion of mass between clay and the liquid in the pores, the equations for the balance of mass for the two constituents are given by
ρ i t + d i v   ( ρ i v i ) = 0 ;    i = s , f
where ( . ) t is the time derivative and d i v is the divergence operator. The balance of linear momentum for the two components are given by
ρ i d i v i d t = d i v   σ i + ρ i b i + i   ;    i = s , f
where b i is the external body force, i is the interaction forces, σ i designates the partial stress tensor and d i ( . ) d t = ( . ) t + [ g r a d   ( . ) ] . v i is the total time derivative. We can define the stress tensor for the mixture (the total stress) as
σ = σ s + σ f
If t s denotes the traction vector of the solid component on the boundary, then
t s = σ s n s
where n s is the unit vector acting on the boundary. In accordance with the basic principles of mixture theory, advocated by Truesdell [26], if we add the two equations for the conservation of mass or the conservation of the linear momentum, we obtain the balance of mass and the linear momentum for the mixture (note that, in this case, i drops out of the equations), due to application of the Newton’s third law
s = f
Finally, the balance of angular momentum indicates that, in the absence of angular momentum supply, the total stress is symmetric. That is [see Rajagopal [27]]
σ = σ T
Since we consider a non-isothermal problem, we do not discuss the balance of energy or the Second Law of thermodynamics.

3. Constitutive Modeling

In this paper, the emphasis is on the modeling of the stress tensor of the solid particles (clay). However, for the closure of the governing equations, we also need constitutive relations for the stress tensor for the fluid, σ f , and the interaction forces i .

3.1. Fluid Component and the Interaction Forces

The following assumptions, as elaborated in Martins-Costa et al. [28] and Rajagopal [27], would lead to the classical Darcy’s equation. The interaction forces designated by i , within the context of mixture theory and many of the multiphase theories, are usually based on generalizing the interaction force for very special cases, such as the Stokes drag on a single spherical particle.
We assume that the frictional (viscous) forces within the fluid can be ignored and as a result the partial stress tensor for the fluid can be given by a Eulerian fluid model:
σ f = p f ( ρ f ) I
where p f is, in general, a function of density, and I is the identity tensor [note that compressive stresses are assumed to be negative in these theories, whereas in, geomechanics-related problems, the opposite convention is used]. We further assume, as is customary in geomechanics problems and basic flows through porous structures (see Oka and Kimoto [29], p. 34):
σ f = φ p I
The interaction force is assumed to be
f = α ( v f v s )
where α is a coefficient that can depend on porosity, viscosity, permeability, etc. This is basically a generalization of the Stokes’ drag on a single spherical particle. In general, the interaction forces could depend on other kinematical quantities such as the relative acceleration, velocity gradients, etc. (see Massoudi [30,31] for a review of this topic). To obtain the Darcy’s equation, we ignore the inertial effects, i.e., we ignore the left-hand side of Equation (18) (when written for the fluid component), and, by using Equations (13) and (23), we obtain the Darcy’s equation:
p f + α w + ρ f b f = 0
where ρ f , as mentioned before, is the density of the fluid. Furthermore, by assuming (see Martins-Costa, et al. [28], Williams [32])
α = μ k
where k represents the specific permeability, which for anisotropic materials is generally a second-order tensor. Equation (26) can be re-written as
w = k μ ( p f ρ f b f )
In soil mechanics literature, Darcy’s equation is sometimes expressed using the concept of hydraulic conductivity ( K ) (see Bear and Bachmat (1990, p. 294)), defined by
K = k ρ f g μ
Here, ρ f g represents the volumetric weight of the fluid, γ f , which is assumed to act in the vertical direction ( b f = [ 0 , γ f , 0 ] T ). Equation (28) is then re-written as:
w = K γ f ( p f ρ f g )
We use this form of the equation in our finite element simulation.

3.2. Solid Component

The partial stress for the solid component can be defined as (Lewis and Schrefler [33]):
σ s = σ ( 1 φ ) p f I
where σ is the effective stress tensor and p f is the pore (fluid) pressure. We can relate the total stress tensor of the mixture σ . [see Equation (19)] to the effective stress tensor σ [by adding Equations (19), (24) and (31)], namely
σ = σ p f I
We assume that the strain in clay can be decomposed into an elastic part and a viscoplastic part. In plasticity theory, we find it more convenient to assume this decomposition applies to for the strain rates [see Davis and Selvadurai [34], p. 97]
ε ˙ = ε ˙ e + ε ˙ v p
We also assume that the elastic part of the strain can be represented by the “small-strain” or the linearized theory of elasticity, where, as customary in soil mechanics, the strain is assumed to depend on the effective stress (Terzaghi [4], pp. 11–15), Schofield and Wroth ([35], p. 9). For an isotropic material, using the index notation, the elastic strain is given by (Matsuoka and Sun [36], p. 37)
ε ˙ i j e = 1 + ν E σ ˙ i j ν E ( σ ˙ k k ) δ i j
Where in accordance with the critical state theory, we assume E = 3 ( 1 2 ν ) ( 1 + e 0 ) p κ is the (modified form of the) Young’s modulus, ν is the Poisson’s ratio, κ is the slope of the unloading and reloading path (see [Matsuoka and Sun [36] (p. 35, Figure 2.8); Desai and Sriwardane [9] (p. 289, Figure 11.7)), p = σ k k 3 is the initial mean pressure, e 0 is the initial void ratio, and δ i j is the Kronecker delta ( δ i j = 1 ( i = j ) ,   0 ( i j ) ). In a more compact form (using the index notation), Equation (34) can be written as
ε ˙ i j e = S i j k l σ ˙ k l
where S i j k l is the fourth-order compliance tensor, related to C i j k l the stiffness tensor. Since we have assumed that the material is isotropic, in short hand notation [recalling Hooke’s law σ ˙ k l = C i j k l ε ˙ k l e ],
S i j k l = { 1 E i = j i , j = 1 , 2 , 3 ν E i j i , j = 1 , 2 , 3 G i = j i , j = 4 , 5 , 6
In Equation (35) σ ˙ k l can be generalized to the case of an elasto-viscoplastic case (see Desai and Sriwardane [9], p. 294, Equation 11.32), i.e.,
σ ˙ k l = D i j k l ( ε ˙ i j ε ˙ i j v p )
where D i j k l is a fourth-order tensor similar to the elastic moduli.
For the viscoplastic modeling of the strain rate, we start with Perzyna, who used the associated flow rule and assumed the plastic potential function coincides with the loading function [20]. To apply this theory to geomaterials, the main challenge is to define the static loading function [12]. By definition, f s , represents the stress state where the strain rate is assumed to be zero. Here, we assume that the viscoplastic part of the strain rate in Equation (33) is based on Perzyna’s approach, where
ε ˙ i j v p = ψ ( F )   f p σ i j
ψ ( F ) = { ψ ( F ) :   F   >   0 0 :   F     0
F = f l f r f r
In the above equations, ψ is the rate sensitivity function, and its functional form can be obtained either experimentally or theoretically; 〈 〉 is the Macaulay’s bracket (in Equation (38), the Macaulay’s bracket ensures that the function inside the bracket only has a value when it is positive, otherwise its magnitude is zero); F is the over-stress function; and f p is a new term, representing the new potential surface (surfaces of the proposed EVP model are obtained by extending the Modified Cam Clay surface [Roscoe and Burland [21]; in our EVP model, we require a total of three surfaces (see Figure 1): the loading surface f l , the reference surface f r , and the potential surface f p ) given by
Potential surface:
f p = p p 2 p c p p p + ( q p M ) 2 = 0
and f l   and   f r are given by the same expression as those in the Perzyna model. They are the dynamic loading function (the potential surface) and the static loading functions, respectively, given by
Reference surface:
f r = p r 2 p c r p r + ( q r M ) 2 = 0
Loading surface:
f l = p l 2 p c l p l + ( q l M ) 2 = 0
where p = σ o c t = σ k k 3 , q = 3 2 τ o c t = [ 3 2 ( σ d ) i j ( σ d ) i j ] 1 2 . Similar expressions for p and q (also known as the deviatoric stress) can also be found in the work by Borja and Kavazanjian [19]. In the principal stress space, σ o c t and τ o c t are the octahedral normal stress and the octahedral shear stress, respectively (see Matsuoka and Sun [36], pp. 29–30). In the above equations, the suffixes r ,   l ,   and   p represent the reference surface, the loading surface, and the potential surface, respectively. The meaning of these surfaces is shown in Figure 1. At any given time, the reference stress state and the reference surface are known from the laboratory test. The current stress state and the potential stress state are related to the reference stress state through the radial mapping rule, which was proposed by Phillips and Sierakowski [37]. In this paper, we used two image parameters and the details are given in Appendix A. It is worth mentioning that Islam and Gnanendran [22] demonstrated the strengths and the limitations of the existing methods using the two-surface approach in the EVP models. They used the associated flow rule for their EVP model, while, in the present paper, we use the non-associated flow rule and the three-surface approach. In the previous paper, the surface shapes are two ellipses, while, in the present study, we assume all the surfaces are given by a single ellipse, which is close to the MCC surface. To formulate the new EVP model, we assume that the “projection center" is in the origin of the stress space, which is identical to the MCC model. However, to define the surfaces, the expression of the slope of the critical state line ( M ) is changed with respect to the b-value ( = σ 2 σ 3 σ 1 σ 3 ) [38] as
M = 6   s i n ϕ b 2 b + 1 3 + ( 2 b 1 ) s i n ϕ
where ϕ is the internal angle of friction at the failure for each b -value test, ranging from 0 (triaxial compression) to 1 (triaxial extension).
We have introduced a new parameter M in order to obtain a more realistic surface shape in the π -plane (see Islam and Gnanendran [22]). It is observed that, in any stress state ( 0 < b - v a l u e 1 ) other than the triaxial compression state ( b - v a l u e = 0 ) , the MCC equivalent surface overestimates the stress. For the sake of completeness, in Figure 2, we compare the EVP model based on the MCC surface with the new modified surface presented in this paper. It is observed that the newly extended MCC surface captures and compares well with the experimental results.
In Figure 3, we can see that, at any arbitrary reference time t ¯ , the soil state is at “A”, where the corresponding void ratio is e ¯ . With time changing from t ¯ to t, due to creep, the soil moves from “A” to “B” where the corresponding void ratio is e . Then, the following expressions are obtained from Borja and Kavazanjian [19]:
e ¯ = e N λ   ln p c l + κ   ln ( p c l p )
e = e N λ   ln p c r + κ   ln ( p c r p )
where λ and κ are the slope of the normal consolidation line (the λ -line) and the unloading-reloading line (the   κ -line), respectively, and e N is the void ratio corresponding to the λ -line when p = 1 kPa at t ¯ .
It should be mentioned that λ , κ and e N are the necessary parameters in the EVP model; these can be obtained either from the oedometer test or the triaxial test. The meaning of these parameters is the same as those in the MCC model. As time increases, the λ -line changes to the λ ˙ -line and e N is transformed to e ˙ N . The λ ˙ -line will generate the new bounding surface. The λ -line and the λ ˙ -line are parallel, as shown by Bjerrum theory [39]. Using Borja and Kavazanjian’s [19] concept and the multisurface theory, t ¯ is the arbitrary time, representing the state of stress prior to the surface evolving. In Equations (40) and (41), p c l and p c r are also known as the creep exclusive preconsolidation pressure and the creep inclusive preconsolidation pressure, respectively (see Islam and Gnanendran [22]). The expression for p c l is similar to the one used in the MCC model:
p c l = p + q 2 p M 2
After rearranging Equation (44), the expression for p c r can be obtained
p c r = e x p ( e N e κ   l n p λ κ )
From Figure 3 for the definitions of λ and κ , we can obtain an expression for p c p :
p c p = p c r λ κ p c l λ κ κ
The detailed derivation of ψ ( F ) and ε ˙ i j v p are presented in Appendix B and Appendix C, respectively. In closing this section, we need to mention that the overstressed EVP models are usually based on the Perzyna’s theory [20] in combination with the critical state soil mechanics theory, e.g., the Modified Cam Clay (MCC) model [21]. In these approaches, the viscous nature of the EVP model is introduced in the theory using a secondary compression index, a creep function and a relaxation function. However, in most cases, to reduce the complexity of the model and to minimize the number of parameters, the EVP models are developed considering the associated flow rule, where it is assumed that the yield surface is identical to the potential surface. To capture the behavior of geomaterials, using the non-associated flow rule is essential [40]. Depending on the application of the Critical States Soil Model (CSSM) in the EVP models using the non-associated flow rule, there are two approaches we can consider: (i) those with critical state [40]; and (ii) those without critical state [41]. The required parameters for the EVP models, satisfying the non-associated flow rule, ranges from 7 parameters [40] to 44 parameters [41]. A summary of the EVP models with the non-associated flow rule for different geomaterials is presented in Table 1.

3.3. Model Parameters

There are six parameters in our model that need to be determined: the consolidation parameters [related to the gradient of the normal consolidation line ( λ ) and the swelling line ( κ ) ], the strength parameter [the slope of the critical state line ( M ) ], the Poisson’s ratio ( ν ) , the state parameter, i.e., the void ratio ( e N ) at the unit mean pressure at any reference time ( t ¯ ) , and the creep parameter, i.e., the secondary compression index ( C α ) . Among these six parameters, five are similar to the MCC model parameters ( λ , κ , M , ν , e N ) . The details of how these parameters can be obtained were given by Roscoe and Burland [21]. The other parameter, C α = Δ e Δ l o g ( t ) , can be obtained from either the oedometer test or the triaxial test. Here, Δ e is the change of the void ratio during the time change Δ l o g ( t ) . In the literature, there are different definitions available for the secondary compression index (see Augustesen et al. [47]). The two most popular definitions are: (i) the void ratio based; and (ii) the strain based. In this paper, we use the first definition. Augustesen et al. [47] (see Figure 8, p. 141), presented a schematic diagram to calculate the secondary compression index. A similar definition can also be found in many soil mechanics textbooks.

3.4. Summary of the Basic Equations and the Assumptions Used in the Code

The equations used in the code and the finite element formulation do not have the same exact forms as those given in Section 3.1 and Section 3.2. In this subsection, we present a summary of the derivation of the equations used in the code and we show how these equations can be obtained from the equations in the earlier sections, if proper assumptions are made. For example, the momentum (equilibrium) equation for the mixture, as used in finite element code, can be obtained if we write Equation (18) for the two components, add them up, assume that the inertial terms can be ignored, use the definition of the total stress tensor [Equation (19)], and the definition for the total density [Equation (11)], then we have
d i v σ + ρ g = 0
where g is the acceleration due to the gravity. If we substitute Equation (32) into Equation (48), and use the convention that a compressive pressure is positive, we have
d i v   ( σ + p f I ) + ρ g = 0
Similarly, if the conservation of mass equations, given by Equation (17), are re-written for the two components, for the cases of constant densities, i.e., when ρ R s = constant   and   ρ R f = constant , they become
φ t + d i v   ( ( 1 φ )   v s ) = 0
φ t + d i v   ( φ v f ) = 0
Adding these two equations, we obtain
d i v   [ φ ( v f v s ) ] + d i v v s = 0
In this equation, the term v r = φ ( v f v s ) is known as the “specific discharge” and is related to the relative velocity w   . We also notice that
d i v   v s = ε v t
[recall that ε v = t r ε was defined as the dilatation (see Equation (6)), and ε v ˙ = ( t r ε ) t = t ( u i x i ) = d i v v s ] . Now, using Equations (53) and (30) in Equation (52), we obtain
d i v   [ K γ f ( p f ρ f g ) ] + ε v t = 0
This expression is also known as the “storage equation” which is the basic equation for the consolidation theory ( see Bear and Bachmat [48], Equation 4.1.46).

4. Finite Element Solutions

The set of equations that need to be solved is: (i) the mass balance equation; (ii) the equilibrium equation; (iii) stress–strain relations; and (iv) strain–displacement relations. The unknown variables are the stresses, the strains, the displacements and the pore pressure. A summary of the equations and the unknown variables is presented in Table 2.
Using the non-associated flow rule EVP model presented in this paper, we can see that there are enough equations to solve for the unknown variables. For any deformable porous media, similar expressions can be found in the work of Bear and Bachmat ([48], Chapter 4).
To solve these equations using finite element, we can use either “the strong form” or “the weak form”. The details of these two techniques can be found in many finite element text books (Zienkiewicz et al. [49]). Even though both approaches provide similar results, considering the “continuity requirements” and the “symmetry of the stiffness matrix”, the weak form solution is better compared to the strong form [49]. In this paper, we use the weak form approach. To obtain the weak form, the Galerkin weighted residual method [50] is introduced. Equations (48) and (54) are used to solve the problem. In the following section, a summary of the weak formulation is described.

Weak Form Formulation

Equilibrium equation
T σ + b = 0     in Ω
Darcy’s equation
w = K γ w ( p f ρ f g   )    in Ω
Continuity equation
T w + m T ε v ˙ = 0     in Ω
Effective stress relation
σ = σ + m p f
In Equations (55)–(58), σ   and   σ are the total stress and the effective stress, respectively, b is the body force, T is the differential operator, T is the divergence operator, is the gradient operator, p f is the pore pressure, w is the superficial velocity, K is the coefficient of permeability, γ w is the unit weight of the fluid, b f is the body force vector for the fluid, and ε ˙ is the volumetric strain rate of soil. The general expression of T , ,   b , b f m , w and   K in the above Equations are
T = [ x 0 0 y 0 z 0 y 0 x z 0 0 0 z 0 y x ]
= [ x , y , z ] T
b = [ b x , b y , b z ] T
b f = [ 0 , γ w , 0 ] T
m = [ 1 , 1 , 1 , 0 , 0 , 0 ] T
w = [ w x , w , w z ] T
K = [ K x 0 0 0 K y 0 0 0 K z ]
Assuming that the soil is subjected to the traction t , the following conditions can be written.
Initial conditions:
u ( t = 0 ) = u 0     in Ω
p f ( t = 0 ) = p f 0     in Ω
Boundary conditions:
for the soil,
σ . n = t     on   Γ t
u = u 0     on   Γ u
for the fluid,
p f = p f 0     on   Γ p f
v . n = 0     on   Γ v
In Equations (66)–(71), u 0 is the initial displacement, p f 0 is the initial fluid pressure, n is the unit normal vector, Γ t , Γ u , Γ p f and Γ v are the respective segment of the boundary, and the subscripts have similar meanings as described above. The surface traction vector is defined a
t = [ t x , t y , t z ] T
t x = σ x n x + τ x y n y + τ x z n z
t y = τ x y n x + σ y n y + τ y z n z
t z = τ x z n x + τ y z n y + σ z n z
In Equations (66)–(71), two types of boundary conditions are used: (i) the essential boundary condition (EBC) or the Dirichlet boundary condition; and (ii) the natural boundary condition (NBC) or the Neumann boundary condition. For the EBC, the primary variables in the domain are known, while, for the NBC, the differential form of the primary variables are prescribed.
Assuming that soil is isotropic, the stress tensor σ , the stain tensor ε and the displacement vector u can be expresses as
σ = [ σ x x , σ y y , σ z z , σ x y , σ y z , σ z x ] T
ε = [ ε x x , ε y y , ε z z , ε x y , ε y z , ε z x ] T
u = [ u x x , u y y , u z z ] T
To obtain the weak formulation for the equilibrium equation, basic necessary steps are: (i) replacing Equation (58) in Equation (55) and multiplying the resulting equation with an arbitrary function, which removes the essential boundary condition or the Dirichlet boundary condition [51]; (ii) integrating the system over the domain; and (iii) applying the Green–Gauss theorem [50] and integrating by parts.
Ω [ B u ] T σ ˙ d Ω + Ω [ B u ] T m p ˙ f d Ω = Ω [ N u ] T b ˙ d Ω + Γ [ N u ] T T ˙ d Γ
Now, σ ˙ (see Equation (37)) can be written as
σ = D ˙ ( B u u ˙ n ) D ε ˙ v p
where, in Equation (37), ϵ ˙ = B u u ˙ n and in Equation (79), p ˙ f = N p p ˙ f n [Zienkiewicz et al. [49]]. The compact form of Equation (79) can be written as
K E u ˙ n + I E p ˙ f n = F E
where
K E = Ω [ B u ] T D B u d Ω
I E = Ω [ B u ] T m N p d Ω
F E = Ω [ N u ] T b ˙ d Ω + Γ [ N u ] T T ˙ d Γ + Ω [ B u ] T D ε ˙ v p d Ω
D = [ ( D e C ) 1 + θ d Δ t n H n ] 1
H n =   ( ε ˙ v p σ ) n
B u = N u
N u = [ N u 1 0 0 ... N u n 0 0 0 N u 1 0 ... 0 N u n 0 0 0 N u 1 ... 0 0 N u n ]
N p = [ N p 1 ... N p n ]
m = [ 1 , 1 , 1 , 0 , 0 , 0 ] T
b = [ b x , b y , b z ] T
= [ x 0 0 0 y 0 0 0 z y x 0 0 z y z 0 x ]
Similar formulation for Equation (57) can be written as
[ I E ] T u ˙ n + H E p f n = Q E
where
[ I E ] T =   Ω [ N p ] T m T B u d Ω
H E = Ω [ B p ] T K γ w B p d Ω
Q E = Ω [ B p ] T K γ w   b f d Ω Γ [ N p ] T w d Γ
B p = N p
= [ x , y , z ] T
Now, combining Equations (81) and (93), we obtain
[ K E I E [ I E ] T 0 ] [ u ˙ n p ˙ f n ] + [ 0 0 0 H E ] [ u n p f n ] = [ F E Q E ]
Equation (99) represent the element matrix. To obtain the global matrix, Equation (99) needs to be summed over a number of elements as follows
K G = n = 1 n K E
I G = n = 1 n I E
[ I G ] T = n = 1 n [ I E ] T
H G = n = 1 n H E
F G = n = 1 n F E
Q G = n = 1 n Q E
In the time interval   Δ t n =   t n + 1 t n , the viscoplastic strain rate increment Δ ε n v p can be defined as [7,8,52]
Δ ε n v p = Δ t n [ ( 1 θ d ) ε ˙ n v p + ε ˙ n + 1 v p ]
where θ d = 0 , denotes the “fully explicit” Euler time integration scheme (or the forward difference method) and θ d = 1 indicates the “fully implicit” Euler time integration scheme (or the backward difference method). For θ d 0.5 , the integration scheme is unconditionally stable and, when θ d < 0.5 , the integration process is conditionally stable. The integration scheme for θ d = 0.5 is the “implicit trapezoidal” scheme or the Crank–Nicolson rule.
The viscoplastic strain rate increment ε ˙ n + 1 v p at time step (n +1) in Equation (106) can also be written as follows [8]
ε ˙ n + 1 v p = ε ˙ n v p + H n Δ σ n
where Δ σ n is the stress change at any time interval Δ t n =   t n + 1 t n and H n (given in Appendix D) is the derivative of the viscoplastic strain-rate with respect to time at the nth time step.
Combining Equations (106) and (107), the above can be written as follows [8]
Δ ε n v p = Δ t n ε ˙ n v p + C n Δ σ n
where,
C n = θ d H n Δ t n
Δ σ n = D ( B n Δ u n ( Δ t n ε ˙ n v p + C n Δ σ n ) )

5. Model Verification and Discussion

We investigated the performance of our EVP model in a variety of experimental applications. These applications included the consolidated undrained triaxial compression test and the extension test, the consolidated drained compression test, the over-consolidation ratio test, the confining pressure effect, the strain rate test, the creep test, and the relaxation test. In this respect, we considered both natural clay and laboratory obtained reconstituted clay. We specifically investigated the Osaka clay [53], the San Francisco Bay Mud (SFBM) clay [54], the Kaolin clay [55], and the Hong Kong Marine Deposit (HKMD) clay [56]. Clay properties are presented in Table 3 and the triaxial loading configuration is shown in Figure 4.

5.1. Simulation of the Undrained Triaxial Test on the Osaka Clay

Adachi et al. [53] presented undrained triaxial compression test for the natural Osaka clay. This is a moderately sensitive clay (sensitivity, S t = 14.50). The water content, the liquid limit and the plastic limit of this clay are 65.00–72.00%, 69.20–75.10% and 24.50–27.30%, respectively. Its specific gravity is 2.62–2.70. The particle size distribution of this clay consists of the clay fraction of 44.00%, the silt fraction of 49.00% and the sand fraction of 7.00%. The strain rate during the test was 1 × 10 2   % / minute . From the experimental results, it is evident that the Osaka clay shows strain softening behavior, which was captured using the new EVP model presented in this paper.
In Figure 5, the experimental results for the Osaka clay for two different mean pressures (176.4 kPa and 235.2 kPa) are compared with the non-associated flow rule type EVP model, and the EVP model developed by Kutter and Sathialingam [42]. In Figure 5a, we can see that, for the mean pressure of 235.2 kPa, the new EVP model and the Kutter and Sathialingam [42] EVP model predictions were identical up to an axial strain of 1.87%. After this point, their model over-predicted the observed experimental results. After an axial strain of 14%, the over-prediction in their model was 16%, whereas the new EVP model captured the strain softening. On the other hand, at an axial strain of 14%, the under-prediction in the new EVP model was 5.20%. For small strain cases (axial strains of 1–3.75%), an over-prediction was observed in both EVP models. Jiang et al. [57] reported that, at small strains, the over-prediction of the nonlinear responses of clay are expected. Using the hysteretic response equation [58], such issues can be resolved; however, this requires additional model parameters. To make the present EVP model formulation simple, and to limit the number of model parameters, the hysteretic response equation concept was not introduced here. A similar trend was also observed for the mean pressure of 176.4 kPa.
In Figure 5b, the stress paths obtained from the associated flow rule, the non-associated flow rule EVP models, and the Kutter and Sathialingam model are compared with the data for the Osaka clay. It is evident that, after the deviatoric stress reached the peak, the attributed stress paths gradually followed the “narrow region”. This phenomenon indicates that the critical state concept is applicable to the natural soft clay, even at large strains [53]. It was observed that the new non-associated flow rule EVP model could capture the “narrow region”, whereas there was marginal under-prediction in the associated flow rule EVP model.

5.2. Simulation of the Undrained Triaxial Stress Relaxation Test on the San Francisco Bay Mud Clay

Lacerda [54] presented the results of undrained triaxial stress relaxation tests for the undisturbed San Francisco Bay Mud clay. In this paper, only SR-I-5 test data are presented and compared with the EVP model predictions. The sample was isotropically consolidated to a pressure of 78.4 kPa. Basic properties of this clay include the specific gravity = 2.66–2.75; the liquid limit = 88.4–90%; the plastic limit = 35–44%; the plasticity index = 45–55%; and the moisture content = 88–93%. Kaliakin and Dafalias [59] reported that the initial void ratio ( e 0 ) and its corresponding mean pressure ( p 0 ) for the San Francisco Bay Mud are 1.30 and 156.9 kPa, respectively. The details of the undrained triaxial stress relaxation tests on the SFBM clay are presented in Table 4.
In Figure 6, the measured data for the San Francisco Bay Mud clay are compared with the numerical simulations for the associated flow rule and the non-associated flow rule EVP models. In addition, our results are also compared with the EVP models of Kaliakin and Dafalias [59] and Kutter and Sathialingam [42]. It is evident that the non-associated flow rule prediction was close to the laboratory results.

5.3. Simulation of the Consolidated Undrained Triaxial Tests on the Kaolin Clay

Herrmann et al. [55] presented the effect of over-consolidation ratio (OCR) on the Kaolin clay for the undrained triaxial compression test and the extension test. This clay is a mixture of the Kaolin (Snow-Cal 50) and a 5.0% Bentonite by weight, which was prepared in the laboratory. The Kaolin clay properties are: the specific gravity = 2.64, the liquid limit = 47.0%, and the plastic limit = 20.0%. The sample was isotropically consolidated to a pressure of 392.2 kPa; to achieve OCR = 1, 2, 4 and 6, the confining pressure was changed accordingly. For the OCR = 1, the corresponding initial void ratio ( e 0 ) is 0.613. The deviatoric stress versus the axial strain predictions for the consolidated undrained triaxial compression test (OCR = 1, 2, 4 and 6) and the extension test (OCR = 1 and 2) were compared with the experimental results, as shown in Figure 7. For the normally consolidated clay (OCR=1), it was observed that the new EVP model captured the stress–strain responses well before the deviatoric stress reached its peak, but then slightly under-predicted them (at 14% strain, the under-prediction was only 1.05 %). In Figure 7a, it is evident that, when OCR = 2, 4 and 6, before reaching the peak deviatoric stress, the model over-predicted but then exhibited small amount of under-prediction. In the consolidated undrained extension tests, a similar trend was observed. When OCR = 1 and 2, we compared the experimental results with the predicted responses of the pore-water pressure for the consolidated undrained triaxial compression test and the consolidated undrained triaxial extension tests, as shown in Figure 7b. For the normally consolidated clay, predictions for both the compression test and the extension test were satisfactory. For the OCR = 2, in the triaxial compression test, the new model slightly over-predicted before the maximum pore pressure was reached and, in the triaxial extension, the negative pore pressure was well captured with a small degree of under-prediction. A similar pattern was observed for the stress path, as shown in Figure 7c.

5.4. Simulation of the Consolidated Drained Tests on the Hong Kong Marine Deposit Clay

The Hong Kong Marine Deposit clay [see Yin and his co-workers [56]] was used to study the isotropically consolidated drained test. This is a reconstituted soft to very soft illitic silty clay, which has been extensively studied.
The properties of the Hong Kong Marine Deposit (HKMD) clay are: the specific gravity = 2.66, the liquid limit = 60%, the plastic limit = 28%, the plasticity index = 32%, and the moisture content = 51.7%. The predictions of the new model for the normally consolidated drained triaxial test on the HKMD clay are shown in Figure 8a for the mean pressures of 300 kPa and 400 kPa. It is evident that the present model captured the deviatoric stress versus axial strain response very well (Figure 8). The volumetric response and the stress paths are presented in Figure 8b,c, which were also satisfactory.

6. Perturbation Analysis

We performed a perturbed analysis for each model parameter for the Osaka clay and the Kaolin clay. The first one is a natural clay, while the second one is a reconstituted clay. The objectives of such analyses were to quantify the sensitivity of each perturbed parameter to the calibrated reference parameters. The mean value of the error ( E m e a n )   was calculated for each parameter as
E   m e a n ( % ) = 1 N n = 1 N | q r e f k q p e r k q r e f k | × 100
where N is the number of error calculations at different strain values; and q r e f k and q p e r k are the predicted values of the deviatoric stress at k t h   strain using a reference parameter and a perturbed model parameter, respectively. The material parameters of the EVP models were divided into two groups: the Modified Cam Clay (MCC) parameters and the viscous parameters. Here, the parameters were perturbed by the same magnitude so that the influence of the perturbation of the parameter on model’s prediction could be observed.
The parameters in the Osaka clay model were perturbed by 3 % of their reference values, which were from experimental data. Then, for each parameter, E m e a n ( % ) due to the perturbation was calculated for the associated flow rule and the non-associated flow rule EVP models, as shown in Equation (111). The perturbed results for the Osaka clay are shown in Figure 9. It was observed that parameters e N and   λ were most sensitive, while M had medium sensitivity and κ   and C α   were less sensitive. To compare the perturbation influence on the flow rule, identical perturbation was also conducted for the Kaolin clay, as shown in Figure 10. In the perturbation analyses with the same magnitude, the effect of the flow rule was negligible.

7. Concluding Remarks

To describe the time-dependent behavior of clay, a new EVP model is formulated in this paper. The model requires six independent parameters, which are defined in the literature and can be extracted from the simple oedometer tests and the triaxial tests. A generalized nonlinear creep function is also used. The model captured a wide range of results when compared with experimental observations obtained from the drained and the undrained conditions for the triaxial compression test and the extension test, the creep test, the relaxation test, the over consolidation ratio effect test, and the strain rate test. Both undisturbed clay and reconstituted clay were considered in our study. From the comparison of the model predictions, it is evident that reasonably good agreement was obtained.

Author Contributions

Conceptualization, M.M.; Investigation, M.N.I.; Methodology, M.N.I. and M.M.; Supervision, C.T.G.; Validation, M.N.I.; Writing—original draft, M.N.I.; Writing—review & editing, M.M.

Funding

This research received no external funding.

Acknowledgments

The first author was financially supported during his research at the University of New South Wales, Canberra, Australia.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of the Image Parameters

In this paper, we use the “radial mapping” rule [37] as discussed by Hashiguchi and Ueno [60] and Kaliakin and Dafalias [59] with the projection center of the yield locus at the center of the stress space (see Figure 1), which is similar to the Modified Cam Clay model (see also Desai and Siriwardane [9] (p. 288, Figure 11.6); Islam and Gnanendran [22] (Appendix III); Kaliakin and Dafalias [59] (Equation (4b)).
To obtain the image stress of the loading surface ( p , q ) (see also Equation (41)) on the reference surface ( p r , q r ) (see also Equation (40)), we use a parameter ( β 1 ) defined as
( σ r ) i j = β 1 σ i j p r   = β 1 p q r   = β 1 q }
Now, substituting Equation (A1) into Equation (40) and solving it for the quadratic real value, we obtain β 1 as
β 1 = p c r p c l
Similarly, for the potential surface we obtain
( σ p ) i j = β 2 σ i j p p   = β 2 p q p   = β 2 q }
β 2 = p c p p c l
where   p c l ,   p c r   and p c p are defined in Equations (45)–(47), respectively.

Appendix B. Derivation of the Rate Sensitivity Function ψ ( F )

Islam and Gnanendran [22] presented the derivation of the rate sensitivity function, ψ ( F ) for the associated flow rule and also discussed the limitations of these approaches. For the sake of completeness, we provide a summary of the derivation for the non-associated flow rule. From Islam and Gnanendran [22], the viscoplastic strain rate ( ε ˙ v v p ) is given as
ε ˙ v v p = α t ¯ ( 1 + e 0 ) ( p c l p c r ) λ κ α
If we decompose ε ˙ i j v p into the volumetric ( ε ˙ v v p ) and the deviatoric ( ε ˙ q v p ) parts, the expression of the viscoplastic strain rate ( ε ˙ i j v p ) in the conventional triaxial stress space can be shown to be
ε ˙ v v p   = ψ f p p p ε ˙ q v p = ψ f p q p }
If we compare Equations (A5) and (A6), the following expression can be obtained for the one-dimensional compression case considering the associated flow rule [22]:
ψ   =   α 0 t ¯ ( 1 + e 0 ) ( p c l p c r ) λ κ α 1 2 p c p [ 1 ς 1 2 ]
ς = p c p p p = 1 + ( η 0 M ) 2
η 0 = 6 ( λ κ ) 2 9 ( λ κ ) 2 + ( 2 λ M ) 2 4 λ

Appendix C. Derivation of the Viscoplastic Strain Rate ε ˙ i j v p

Recalling Equation (38), the viscoplastic strain rate ( ε ˙ i j v p ) can be written as
ε ˙ i j v p   =   ψ ( F )   f p σ i j
The expression for ψ ( F ) can be obtained from Appendix B. The mathematical form of f p σ i j can be presented in terms of the chain rule:
f p σ i j p = f p p p p p σ i j p + f p q p q p σ i j p + f p M M b b σ i j p
In Equation (A9), the expression of the third term on the right-hand side depends on the definition of M, which is discussed in Section 3. The different terms in Equation (A9) can be expressed as
p p σ i j p = 1 3 δ i j ,   w h e r e   δ i j = { 1     i f         i = j 0       i f     i j
q p σ i j p = { 3 2 q p ( σ i j p p p δ i j )     i f       i = j 3 2 q p ( 2 σ i j p )     i f     i j
f p M = 2 q p 2 M 3
M b = 3 s i n ϕ ( 2 b 1 ) ( b 2 b + 1 ) [ 3 + ( 2 b 1 ) s i n ϕ ] 12 s i n 2 ϕ ( b 2 b + 1 ) [ 3 + ( 2 b 1 ) s i n ϕ ] 2
b σ 11 p = σ 22 p σ 33 p ( σ 11 p σ 33 p ) 2
b σ 22 p = 1 σ 11 p σ 33 p
b σ 33 p = 1 σ 11 p σ 33 p + σ 22 p σ 33 p ( σ 11 p σ 33 p ) 2
In addition, following Lade [61], the third term in Equation (A9), also can be expressed in terms of Lode angle (θ) as follows
M θ = 1 2 k   s i n ( 3 θ ) c o s ( π 6 + 1 3 c o s 1 ( c o s ( 3 θ ) ) k 27 ) s i n ( π 6 + 1 3 c o s 1 ( c o s ( 3 θ ) ) k 27 ) 2 3 c o s ( 3 θ ) 2 k + 81
θ σ 1 p = 1 3 c o s 1 [ 1 2 1 σ 1 p σ 3 p ( 6 b 3 + 6 b 2 + 3 b ) ( b 2 b + 1 ) 3 2 3 4 1 ( b 2 b + 1 ) 5 2 { ( 2 b 3 3 b 2 3 b + 2 ) ( 1 σ 1 p σ 3 p ( 2 b 2 + b ) ) } ]
θ σ 2 p = 1 3 c o s 1 [ 1 2 1 σ 1 p σ 3 p ( 6 b 2 6 b 3 ) ( b 2 b + 1 ) 3 2 3 4 1 ( b 2 b + 1 ) 5 2 { ( 2 b 3 3 b 2 3 b + 2 ) ( 1 σ 1 p σ 3 p ( 2 b 1 ) ) } ]
θ σ 3 p = 1 3 c o s 1 [ 1 2 1 σ 1 p σ 3 p ( 6 b 3 12 b 2 + 3 b + 3 ) ( b 2 b + 1 ) 3 2 3 4 1 ( b 2 b + 1 ) 5 2 { ( 2 b 3 3 b 2 3 b + 2 ) ( 1 σ 1 p σ 3 p ( 2 b 2 3 b + 1 ) ) } ]

Appendix D. Derivation of the Gradient Matrix H n

Owen and Hinton [8] presented the definition of the gradient matrix as
H n =   ( ε ˙ v p σ ) n
where n is the number of dimensions. For two-dimensional cases, the total number of elements are sixteen, which can be expressed as
[ H ] = [ ε ˙ 11 v p σ 11 ε ˙ 11 v p σ 22 ε ˙ 11 v p σ 33 ε ˙ 11 v p σ 12 ε ˙ 22 v p σ 11 ε ˙ 22 v p σ 22 ε ˙ 22 v p σ 33 ε ˙ 22 v p σ 12 ε ˙ 33 v p σ 11 ε ˙ 33 v p σ 22 ε ˙ 33 v p σ 33 ε ˙ 33 v p σ 12 ε ˙ 12 v p σ 11 ε ˙ 12 v p σ 22 ε ˙ 12 v p σ 33 ε ˙ 12 v p σ 12 ]
Applying the chain rule to the derivatives ( ε ˙ i j v p σ i j ) of the strain rate with respect to the stress components, we have
( ε ˙ i j v p σ i j = ε ˙ i j v p p p σ i j + ε ˙ i j v p q q σ i j )
Implementing Equation (A12) into Equation (A11), we obtain the complete form of the gradient matrix for two-dimensional case:
[ H ] = [ ε ˙ 11 v p p p σ 11 +   ε ˙ 11 v p q q σ 11 ε ˙ 11 v p p p σ 22 +   ε ˙ 11 v p q q σ 22 ε ˙ 11 v p p p σ 33 +   ε ˙ 11 v p q q σ 33 ε ˙ 11 v p p p σ 12 +   ε ˙ 11 v p q q σ 12 ε ˙ 22 v p p p σ 11 +   ε ˙ 22 v p q q σ 11 ε ˙ 22 v p p p σ 22 +   ε ˙ 22 v p q q σ 22 ε ˙ 22 v p p p σ 33 +   ε ˙ 22 v p q q σ 33 ε ˙ 22 v p p p σ 12 +   ε ˙ 22 v p q q σ 12 ε ˙ 33 v p p p σ 11 +   ε ˙ 33 v p q q σ 11 ε ˙ 33 v p p p σ 22 +   ε ˙ 33 v p q q σ 22 ε ˙ 33 v p p p σ 33 +   ε ˙ 33 v p q q σ 33 ε ˙ 33 v p p p σ 12 +   ε ˙ 33 v p q q σ 12 ε ˙ 12 v p p p σ 11 +   ε ˙ 12 v p q q σ 11 ε ˙ 12 v p p p σ 22 +   ε ˙ 12 v p q q σ 22 ε ˙ 12 v p p p σ 33 +   ε ˙ 12 v p q q σ 33 ε ˙ 12 v p p p σ 12 +   ε ˙ 12 v p q q σ 12 ]

References

  1. Saito, M.; Uezawa, H. Failure of soil due to creep. In Proceedings of the 5th International Conference on Soil Mechanics and Foundation Engineering (ICSMFE), Paris, France, 17–22 July 1961; Volume 1, pp. 315–318. [Google Scholar]
  2. Muller-Salzburg, L. The Rock Slide in the Vaiont Valley; Springer: Wien, Germany, 1964. [Google Scholar]
  3. Rowe, P.W. The stress-dilatancy relation for static equilibrium of an assembly of particles in contact. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1962, 269, 500–527. [Google Scholar] [CrossRef]
  4. Terzaghi, K. Theoretical Soil Mechanics; Wiley: New York, NY, USA, 1943. [Google Scholar]
  5. Barden, L. Consolidation of Clay with Non-linear Viscosity. Géotechnique 1965, 15, 345–362. [Google Scholar] [CrossRef]
  6. Alonso, E.E.; Gens, A.; Lloret, A. Precompression design for secondary settlement reduction. Geotechnique 2000, 50, 645–656. [Google Scholar] [CrossRef]
  7. Potts, D.M.; Zdravkovic, L. Finite Element Analysis in Geotechnical Engineering: Theory; Thomas Telford: London, UK, 1999. [Google Scholar]
  8. Owen, D.R.J.; Hinton, E. Finite Elements in Plasticity: Theory and Practice; Pineridge Press Limited: Swansea, UK, 1980. [Google Scholar]
  9. Desai, C.; Siriwardane, H.J. Constitutive Laws For Engineering Materials; Prentice-Hall: Upper Saddle River, NJ, USA, 1984. [Google Scholar]
  10. Chaboche, J.L. A review of some plasticity and viscoplasticity constitutive theories. Int. J. Plast. 2008, 24, 1642–1693. [Google Scholar] [CrossRef]
  11. Tatsuoka, F.; Bendetto, H.D.; Enomoto, T.; Kawabe, S.; Kongkitkul, W. Various viscosity types of geomaterials in shear and their mathematical expression. Soils Found. 2008, 48, 41–60. [Google Scholar] [CrossRef]
  12. Liingaard, M.; Augustesen, A.; Lade, P. Characterization of Models for Time-Dependent Behavior of Soils. Int. J. Geomech. 2004, 4, 157–177. [Google Scholar] [CrossRef]
  13. Singh, A.; Mitchell, J.K. General Stress-strain-time Function for Soils. J. Soil Mech. Found. Div ASCE 1968, 94, 21–46. [Google Scholar]
  14. Feda, J. Creep of Soils: And Related Phenomena (Developments in Geotechnical Engineering), 2nd ed.; Elsevier Science: Amsterdam, The Netherlands, 1992. [Google Scholar]
  15. Adachi, T.; Okano, M. A constitutive equation for normally consolidated clay. Soils Found. 1974, 14, 55–73. [Google Scholar] [CrossRef]
  16. Sekiguchi, H. Macrometric Approaches-Static-Intricsically Time Dependent Constitutive Laws of Soils. In Proceedings of the 11th ICFMFE, San Francisco, CA, USA, 29–31 March 2019; pp. 66–98. [Google Scholar]
  17. Sekiguchi, H. Rheological Characteristics of Clays. In Proceedings of the 9th International Conference on Soil Mechanics and Foundation Engineering, Tokyo, Japan, 10–15 July 1977; pp. 289–292. [Google Scholar]
  18. Dafalias, Y.F.; Herrmann, L.R. Bounding surface formulation of soil plasticity. In Soil Mechanics-Transient and Cyclic Loads; Pande, G.N., Zienkiewicz, O.C., Eds.; Wiley: Chichester, UK, 1982; pp. 253–282. [Google Scholar]
  19. Borja, R.I.; Kavazanjian, J.E. A constitutive model for the stress-strain&-time behaviour of ‘wet’ clays. Geotechnique 1985, 35, 283–298. [Google Scholar]
  20. Perzyna, P. Constitutive equations for rate-sensitive plastic materials. Q. Appl. Math. 1963, 20, 321–331. [Google Scholar] [CrossRef]
  21. Roscoe, K.H.; Burland, J.B. On the generalized stress-strain behavior of wet clay. In Engineering Plasticity; Heyman, J., Leckie, F.A., Eds.; Cambridge University Press: Cambridge, UK, 1968; pp. 535–609. [Google Scholar]
  22. Islam, M.N.; Gnanendran, C.T. Elastic-Viscoplastic Model for Clays: Development, Validation, and Application. J. Eng. Mech. 2017, 143. [Google Scholar] [CrossRef]
  23. Carter, J.P.; Balaam, N.P. AFENA User’s Manual. Version 5.0; Center for Geotechnical Research, University of Sydney: Sydney, Australia, 1995. [Google Scholar]
  24. Atkin, R.J.; Craine, R.E. Continuum theories of mixtures: Basic theory and historical development. Q. J. Mech. Appl. Math. 1976, 29, 209–244. [Google Scholar] [CrossRef]
  25. Bowen, R.M. Theory of Mixtures, Part I. In Continuum Physics III; Eringen, A., Ed.; Academic Press: Cambridge, MA, USA, 1976; pp. 2–127. [Google Scholar]
  26. Truesdell, C. Rational Thermodynamics, 2nd ed.; Springer: New York, NY, USA, 1984. [Google Scholar]
  27. Rajagopal, K.R. On a hierarchy of approximate models for flows of incompressible fluids through porous solids. Math. Models Methods Appl. Sci. 2007, 17, 215–252. [Google Scholar] [CrossRef]
  28. Martins-Costa, M.L.; Sampaio, R.; da Gama, R.M.S. Modelling and simulation of energy transfer in a saturated flow through a porous medium. Appl. Math. Model. 1992, 16, 589–597. [Google Scholar] [CrossRef]
  29. Oka, F.; Kimoto, S. Computational Modeling of Multiphase Geomaterials; CRC Press: London, UK, 2013. [Google Scholar]
  30. Massoudi, M. On the importance of material frame-indifference and lift forces in multiphase flows. Chem. Eng. Sci. 2002, 57, 3687–3701. [Google Scholar] [CrossRef]
  31. Massoudi, M. Constitutive relations for the interaction force in multicomponent particulate flows. Int. J. Non-Linear Mech. 2003, 38, 313–336. [Google Scholar] [CrossRef]
  32. Williams, W.O. Constitutive equations for flow of an incompressible viscous fluid through a porous medium. Q. Appl. Math. 1978, 36, 255–267. [Google Scholar] [CrossRef]
  33. Lewis, R.W.; Schrefler, B.A. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, 2nd ed.; John Wiley & Sons: New York, NY, USA, 1998. [Google Scholar]
  34. Davis, R.O.; Selvadurai, A.P.S. Plasticity and Geomechanics; Cambridge University Press: New York, NY, USA, 2002. [Google Scholar]
  35. Schofield, A.N.; Wroth, P. Critical State Soil Mechanics; McGrawHill: London, UK, 1968. [Google Scholar]
  36. Matsuoka, H.; Sun, D.A. The SMP Concept-Based 3D Constitutive Models for Geomaterials; Taylor & Francis: New York, NY, USA, 2006. [Google Scholar]
  37. Phillips, A.; Sierakowski, R.L. On the concept of the yield surface. Acta Mech. 1965, 1, 29–35. [Google Scholar] [CrossRef]
  38. Habib, P. Influence of the variation of the average principal stress upon the shearing strength of soils. In Proceedings of the 3rd International Conference on Soil Mechanics and Foundation Engineering, Zurich, Switzerland, 16–27 August 1953; pp. 131–136. [Google Scholar]
  39. Bjerrum, L. Engineering geology of Norwegian normally consolidated marine clays as related to settlements of buildings. Geotechnique 1967, 17, 81–118. [Google Scholar] [CrossRef]
  40. Zienkiewicz, O.C.; Humpheson, C.; Lewis, R.W. Associated and non-associated visco-plasticity and plasticity in soil mechanics. Geotechnique 1975, 25, 671–689. [Google Scholar] [CrossRef]
  41. Maranini, E.; Yamaguchi, T. A non-associated viscoplastic model for the behaviour of granite in triaxial compression. Mech. Mater. 2001, 33, 283–293. [Google Scholar] [CrossRef]
  42. Kutter, B.L.; Sathialingam, N. Elastic-viscoplastic modelling of the rate-dependent behaviour of clays. Géotechnique 1992, 42, 427–441. [Google Scholar] [CrossRef]
  43. Hickman, R.J.; Gutierrez, M.S. Formulation of a three-dimensional rate-dependent constitutive model for chalk and porous rocks. Int. J. Numer. Anal. Methods Geomech. 2007, 31, 583–605. [Google Scholar] [CrossRef]
  44. Cristescu, N.D. Nonassociated elastic/viscoplastic constitutive equations for sand. Int. J. Plast. 1991, 7, 41–64. [Google Scholar] [CrossRef]
  45. Florea, D. Nonassociated elastic/viscoplastic model for bituminous concrete. Int. J. Eng. Sci. 1994, 32, 87–93. [Google Scholar] [CrossRef]
  46. Jin, J.; Cristescu, N.D. An elastic/viscoplastic model for transient creep of rock salt. Int. J. Plast. 1998, 14, 85–107. [Google Scholar] [CrossRef]
  47. Augustesen, A.; Liingaard, M.; Lade, P. Evaluation of Time-Dependent Behavior of Soils. Int. J. Geomech. 2004, 4, 137–156. [Google Scholar] [CrossRef]
  48. Bear, J.; Bachmat, Y. Introduction to Modeling of Transport Phenomena in Porous Media; Kluwer Academic Publishers: London, UK, 1990; Volume 4. [Google Scholar]
  49. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals, 6th ed.; Elsevier: New York, NY, USA, 2005. [Google Scholar]
  50. Galerkin, B.G. Series solution of some problems of elastic equilibrium of rods and plates (in Russian). Vestn. Inzh. Tech. 1915, 19, 897–908. [Google Scholar]
  51. Koutromanos, I. Fundamentals of Finite Element Analysis: Linear Finite Element Analysis; John Wiley and Sons: Noida, India, 2018. [Google Scholar]
  52. Kanchi, M.B.; Zienkiewicz, O.C.; Owen, D.R.J. The visco-plastic approach to problems of plasticity and creep involving geometric non-linear effects. Int. J. Numer. Methods Eng. 1978, 12, 169–181. [Google Scholar] [CrossRef]
  53. Adachi, T.; Oka, F.; Hirata, T.; Hashimoto, T.; Nagaya, J.; Mimura, M.; Pradhan, T.B.S. Stress–strain behavior and yielding characteristics of Eastern Osaka clay. Soil Found. 1995, 35, 1–13. [Google Scholar] [CrossRef]
  54. Lacerda, W.A. Stress-Relaxation and Creep Effects on Soil Deformation. Ph.D. Thesis, University of California, Berkeley, CA, USA, 1976. [Google Scholar]
  55. Herrmann, L.R.; Shen, C.K.; Jafroudi, S.; DeNatale, J.S.; Dafalias, Y.F. A Verification Study for the Bounding Surface Plasticity Model for Cohesive Soils; Final Report to the Civil Engineering Laboratory; Naval Construction Battalion Center: Port Hueneme, CA, USA, 1982. [Google Scholar]
  56. Yin, J.H.; Zhu, J.G. Measured and predicted time-dependent stress-strain behaviour of Hong Kong marine deposits. Can. Geotech. J. 1999, 36, 760–766. [Google Scholar] [CrossRef]
  57. Jiang, J.; Ling, H.I.; Kaliakin, V.N. An Associative and Non-Associative Anisotropic Bounding Surface Model for Clay. J. Appl. Mech. 2012, 79, 031010. [Google Scholar] [CrossRef]
  58. Whittle, A.; Kavvadas, M. Formulation of MIT-E3 Constitutive Model for Overconsolidated Clays. J. Geotech. Eng. 1994, 120, 173–198. [Google Scholar] [CrossRef]
  59. Kaliakin, V.N.; Dafalias, Y.F. Verification of the Elastoplastic-Viscopalstic Bounding Surface Model for Cohesive Soils. Soils Found. 1990, 30, 25–36. [Google Scholar] [CrossRef]
  60. Hashiguchi, K.; Ueno, M. Elastoplastic constitutive laws of granular materials. In Proceedings of the 9th International Conference on Soil Mechanics and Foundation Engineering (ICFSME), Session 9, Tokyo, Japan, 10–15 July 1977; pp. 73–82. [Google Scholar]
  61. Lade, P.V. Triaxial Testing of Soils; Chapter 2; John Wiley & Sons: Chichester, UK, 2016; p. 88. [Google Scholar]
Figure 1. Illustration of the reference surface, the loading surface, and the potential surface in the p - q plane.
Figure 1. Illustration of the reference surface, the loading surface, and the potential surface in the p - q plane.
Geosciences 09 00145 g001
Figure 2. Comparison of the MCC surface and the extended surface in the π -plane.
Figure 2. Comparison of the MCC surface and the extended surface in the π -plane.
Geosciences 09 00145 g002
Figure 3. Relative locations of p c l , p c r and p c p .
Figure 3. Relative locations of p c l , p c r and p c p .
Geosciences 09 00145 g003
Figure 4. The stress condition in triaxial compression/extension tests (a) cylindrical sample, (b) block sample, (c) axisymmetric/biaxial representation of solid and fluids in the sample.
Figure 4. The stress condition in triaxial compression/extension tests (a) cylindrical sample, (b) block sample, (c) axisymmetric/biaxial representation of solid and fluids in the sample.
Geosciences 09 00145 g004
Figure 5. Comparison of the measured and the predicted consolidated undrained triaxial compression tests results for the Osaka clay: (a) the deviatoric stress vs. axial strain; and (b) the stress path.
Figure 5. Comparison of the measured and the predicted consolidated undrained triaxial compression tests results for the Osaka clay: (a) the deviatoric stress vs. axial strain; and (b) the stress path.
Geosciences 09 00145 g005
Figure 6. Comparison of the measured and the predicted undrained triaxial tests for the stage-changed axial strain rate combined with stress relaxation on the SFBM clay: (a) the deviatoric stress vs. axial strain; and (b) the pore pressure vs. axial strain.
Figure 6. Comparison of the measured and the predicted undrained triaxial tests for the stage-changed axial strain rate combined with stress relaxation on the SFBM clay: (a) the deviatoric stress vs. axial strain; and (b) the pore pressure vs. axial strain.
Geosciences 09 00145 g006
Figure 7. Comparison of the measured and the predicted consolidated undrained triaxial compression and the extension tests on the Kaolin clay: (a) the deviatoric stress vs. the axial strain; (b) the pore-water pressure vs. the axial strain; and (c) the stress path.
Figure 7. Comparison of the measured and the predicted consolidated undrained triaxial compression and the extension tests on the Kaolin clay: (a) the deviatoric stress vs. the axial strain; (b) the pore-water pressure vs. the axial strain; and (c) the stress path.
Geosciences 09 00145 g007
Figure 8. Comparison of the measured and the predicted consolidated drained triaxial compression tests on the HKMD clay: (a) the deviatoric stress vs. the axial strain; (b) the volumetric strain vs. the axial strain; and (c) the stress path.
Figure 8. Comparison of the measured and the predicted consolidated drained triaxial compression tests on the HKMD clay: (a) the deviatoric stress vs. the axial strain; (b) the volumetric strain vs. the axial strain; and (c) the stress path.
Geosciences 09 00145 g008
Figure 9. Sensitivity of perturbation for Osaka clay.
Figure 9. Sensitivity of perturbation for Osaka clay.
Geosciences 09 00145 g009
Figure 10. Sensitivity of perturbation for the Kaolin clay.
Figure 10. Sensitivity of perturbation for the Kaolin clay.
Geosciences 09 00145 g010
Table 1. Summary of the EVP models using the non-associated flow rule for different geomaterials.
Table 1. Summary of the EVP models using the non-associated flow rule for different geomaterials.
ModelMaterialParameterCSSMSurface Shape
Model in this paperClay6YesNC
Zienkiewicz et al. [40]Clay9YesNC
Kutter and Sathialingam [42]Clay7YesC
Hickman and Gutierrez [43]Chalk10YesNC
Cristescu [44]Sand15NoC
Florea [45]Concrete15NoC
Jin and Cristescu [46]Rock32NoC
Maranini and Yamaguchi [41]Granite44NoC
Note: C, Circular; NC, Non-circular; CSSM, Critical state soil mechanics.
Table 2. Balance equations, constitutive relations and the unknown variables.
Table 2. Balance equations, constitutive relations and the unknown variables.
Equations Unknown Variables
NameReferenceNumberVariablesNumber
Mass balanceEquation (54)1 p f 1
Equilibrium Equation (49)3 σ 6
Stress–strain Equation (37)6 ε 6
Strain–displacementEquation (7)6 u 3
Total equation 16Total unknowns16
Table 3. Model parameters.
Table 3. Model parameters.
ParametersOsaka Clay [53]SFBM Clay [54]Kaolin Clay [55]HKMD Clay [56]
λ 0.360.370.150.20
κ 0.0470.0540.0180.045
M b = 0 1.281.401.251.26
M b = 1 0.950.89
ν 0.3G = 23520 kPa0.300.30
e N 3.543.171.512.18
C α 0.03270.0530.0140.106
M b = 0 and M b = 1 are for the triaxial compression and the extension test respectively.
Table 4. Stress relaxation test on the SFBM clay.
Table 4. Stress relaxation test on the SFBM clay.
Phase12345678
TestShearRel.ShearRel.ShearRel.ShearRel.
ε ˙ a 1.501.500.016200.000810
Time0.2530701.281320101.24270016798370
ε a ( % ) 0-0.380.380.38-2.32.32.3-3.943.943.94-5.35.3
Rel. = Relaxation, ε ˙ a = ( % /   min ) , Time = minutes.

Share and Cite

MDPI and ACS Style

Islam, M.N.; Gnanendran, C.T.; Massoudi, M. Finite Element Simulations of an Elasto-Viscoplastic Model for Clay. Geosciences 2019, 9, 145. https://doi.org/10.3390/geosciences9030145

AMA Style

Islam MN, Gnanendran CT, Massoudi M. Finite Element Simulations of an Elasto-Viscoplastic Model for Clay. Geosciences. 2019; 9(3):145. https://doi.org/10.3390/geosciences9030145

Chicago/Turabian Style

Islam, Mohammad N., Carthigesu T. Gnanendran, and Mehrdad Massoudi. 2019. "Finite Element Simulations of an Elasto-Viscoplastic Model for Clay" Geosciences 9, no. 3: 145. https://doi.org/10.3390/geosciences9030145

APA Style

Islam, M. N., Gnanendran, C. T., & Massoudi, M. (2019). Finite Element Simulations of an Elasto-Viscoplastic Model for Clay. Geosciences, 9(3), 145. https://doi.org/10.3390/geosciences9030145

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop