Modification of a Constitutive Model in the Framework of a Multilaminate Method for Post-Liquefaction Sand

Abstract

The estimation of the value and direction of post-liquefaction deformations is one of the most challenging issues in the modelling of liquefaction soil, due to the inherent and induced anisotropy. It is very important in the science of soil-constitutive models to present a simple and comprehensive model for the prediction of fabric anisotropy effects in pre- and post-liquefaction behaviour in granular soil. In the framework of the multilaminate method, 17 planes with pre-determined directions are defined, instead of defining all occurrences depending on the direction in three planes perpendicular to each other in a Cartesian coordinate system. As a result, calculation accuracy is increased in the point due to the effectiveness of the behaviours in different directions. In the present study, after modifying an advanced model by removing constants related to the fabric effect and using lower constants, the precision of model performance after the removal of constants was studied and compared with experimental results in different monotonic, cyclic, drained, and undrained loading conditions. After this, the formation of stress and strain in 17 planes was evaluated in terms of pre- and post-liquefaction, with monotonic and cyclic loadings. The study of the curves shows induced anisotropy in different directions of sandy soil and thus proves the capability of the model in this regard.

Keywords:
multi-directional; granular materials; anisotropy; undrained; liquefaction; cyclic

1 INTRODUCTION

Saturation soil without cohesion undergoes flow failure with huge displacements during rapid static loading or earthquakes. Therefore, it is very important to study the post-liquefaction conditions of sandy soil, due to their significant effects on related structures. Some researchers have evaluated post-liquefaction displacements by using a set of experimental and centrifuge tests in different conditions (Ishikawa et al., 2015; Wang et al., 2015; Basari and Ozden, 2013). Wang et al. (2015) studied post-liquefaction soil in cyclic loading under different post-liquefaction reconsolidations. Their study was conducted on silt soil. It has been observed that if post-liquefaction soil has undergone monotonic loading, soil shear strength and stiffness is increases due to increasing post-liquefaction reconsolidation percent. In addition, it has been observed that the critical state line is different in pre- and post-liquefaction stages. Ishikawa et al. (2015) studied post-liquefaction progressive failure of shallow foundations in centrifuge tests. An important point in induced deformation resulting from liquefaction is its slow movement (several minutes) compared to the shaking duration of an earthquake (several tens of milliseconds). They proved that post-liquefaction progressive failure in shallow foundations results from the spread of the local shear boundaries of the liquefaction layer. Settlement and tilting occur in complete liquefaction, even after the earthquake. Shamoto et al. (1998) divided deformations resulting from liquefaction into two categories: volumetric deformations (which caused settlement) and shear deformations (which caused lateral spreading). They presented a method to simultaneously predict these two displacements in terms of their effect on each other.

A few researchers have studied issues related to post-liquefaction conditions of sand particles using numerical models (Elgmal et al., 2002; Zang and Wang, 2012; Sadrnejad, 2007). Elgmal et al. (2002) introduced a constitutive model to predict liquefaction and numerically studied the effect of the frequency content of input excitation on the post-liquefaction shear deformation. They proved that dominant excitation frequency had the most significant effect on post-liquefaction soil response. Low excitation frequency creates many lateral deformations. Post-liquefaction shear deformation is one of the challenging subjects in liquefaction soil modelling. Zang and Wang (2012) presented a numerical algorithm and its application in a theoretical framework to predict post-liquefaction deformations of saturation sand under cyclic and undrained loadings. They considered the post-liquefaction mechanism by decomposition of volumetric strain into three components with a certain physical background. The interaction between these components controls the post-liquefaction deformations and determines changes in these three physical states in the liquefaction process. This assumption determines the complex issue of transfer of small pre-liquefaction displacements to large post-liquefaction displacements.

Sand particles are formed under different environmental conditions during their lifetime. Bedding plane gradient leads to the formation of weak and strong planes in different directions. The initial anisotropy causes different behaviour of sand in different loading conditions. Another type of anisotropy seen in granular materials is under loading effect and the production of plastic strain in different directions. This type of anisotropy is known as induced anisotropy. There are important effects on the behaviour of sandy soil in different conditions. The three important features of granular soil are liquefaction, dilation, and critical state that can be produced during one loading; these are all affected by induced anisotropy. The post-liquefaction behaviour of sand particles and the direction of their deformations are affected by initial and induced anisotropy, due to the effect of initial anisotropy and the history of stress and strain.

Fabric anisotropic effects were studied on sandy soil responses by a series of experimental studies (Yu et al., 2013; Guo and Zhao, 2013; Louis et al., 2009). Some tried to model important properties of sand particles, such as transfer, liquefaction, and critical state, by determining the relationship between microscopic and macroscopic structures (Guo and Zhao, 2013; Yang et al., 2008; Yang et al., 2013; Hicher and Chang, 2007; Kruyt and Rothenburg, 2016; Khalili and Mahboubi, 2014;Voiadjis et al., 1995; Kruyt and Rothenburg, 2014; Guo, 2014; Li and Yu, 2014; Yan, 2009; Chang and Hicher, 2005). Other studies have been conducted to predict the fabric anisotropic behaviour using different parameters and placing them in equations of constitutive model of soil (Yu et al., 2013; Yang et al., 2008; Gao et al., 2014; Dafalias and Manzari, 2004; Lashkari, 2009; Zhao and Guo, 2013; Dafalias et al., 2004). In addition, a few studies have been conducted to consider the effect of induced anisotropy by using constitutive models (Hareb and Doanh, 2012; Tang Tron Tran et al., 2014; Ye et al., 2012). Although various models have been presented to predict sandy soil behaviour in different loading conditions (Zang and Wang, 2012; Yang et al., 2008; Dafalias and Manzari, 2004; Lashkari, 2009; Liu et al., 2016; Wang and Xie, 2014; Li, 2002; Wang et al., 2014), such models either have been used in restrictive anisotropic conditions and limited loading conditions or they have used a high number of material constants. Most models are not able to determine the parameters dependent on direction, such as fabric, by using stress and strain invariants. Therefore, the multilaminate constitutive model-also known as micro-plane in some sources (Bazant et al., 1996; Caner and Bazant, 2000)-is a suitable mechanism for the simultaneous prediction of complete anisotropy in sandy soil in terms of material behaviour in different planes. This model has been evaluated for different materials, such as soil (Sadrnejad et al., 2009; Sadrnejad, 2009; Sadrnejad and Karimpour, 2011), and for the estimation of damage models in concrete (Ghadrdan and Sadrnejad, 2015; Sadrnejad and Labibzadeh, 2006). Some researchers considered 13 planes in different directions of a sphere with a radius of one; they modelled effects of anisotropy for different models (Sadrnejad, 2009; Sadrnejad and Karimpour, 2011; Sadrnejad and Labibzadeh, 2006). The recent model has been problematic due to the number of planes, preparation of simultaneous conditions, strain consistency, and equilibrium (Ghadrdan and Sadrnejad, 2015; Sadrnejad and Labibzadeh, 2006).

In the present research, the weight coefficients and their directions for 17 planes defined in different directions have been calculated using multilaminate theory; the effects of 17 planes are transferred to points in the framework of the numerical integration method. For constitutive modelling in planes, numerically advanced model, modified with the reduction in number of soil constants, have been used based on limit state and hypo-plasticity theory (Manzari and Dafalias, 1997; Dafalias,1986; Wang et al., 1990). The model should be able to predict the behaviour of sand particles in drained, undrained, cyclic, and monotonic conditions. In the framework of the multilaminate method, the governing equations the plane are presented briefly and the model constants are obtained in planes with calibration. The model capability is proved by studying and comparing the removal of constants from the initial model (Dafalias and Manzari, 2004) using the multilaminate framework and comparing the numerical results of the model to experimental results in different loadings, confining pressure, and void ratio. In the following section, important effects of induced anisotropy are studied in 17 planes in different loading conditions to determine the post-liquefaction conditions of sand particles. The effects of induced anisotropy are studied by plotting different curves in planes. Active, inactive, and progressive planes in failure are determined in triaxial standard tests. Then, the model capability is evaluated in terms of the prediction of behaviour in different directions and the determination of movement path of granular soil after liquefaction.

2 MULTILAMINATE MODEL

The multilaminate theory (Sadrnejad, 2011) is based on the determination of numerical relationship between micro-scale behaviours and engineering mechanical properties (macro scale behaviour) in form of constitutive equation. In other words, material properties are obtained by features of each constitutive element and stress-strain behaviour of material is obtained by studying micro scale behaviour.

When a multi-dimensional part undergoes small shear stress, it carries elastic shear deformation. Multi-dimensional parts start moving in the direction of boundary planes-known as sliding planes-due to the increase in stress and reaching to a certain amount of stress. Shear stress that creates higher deformation is increased by increasing deformation. The total shear deformation is the sum of elastic shear deformation in multi-dimensional parts and plastic shear deformation resulting from the sliding of adjacent parts. When stress is decreased, the elasticity component returns to the starting point. Then, the multi-dimensional parts start sliding inversely due to the reduction in stress, and reaches a certain value. The shear stress required for sliding depends on the normal stress. Sliding occurs when the stress state crosses the yield limit. However, sliding occurs only in sliding planes in the suggested directions (schematic view in Figure 1). On the basis of this, the higher the number of pre-determined planes, the closer to the reality will be the sliding, opening, and closing of the planes.

Figure 1:
(a) Presentation of real accumulation of particles (b) 2D presentation of accumulation of artificial multi-dimensional parts.

In the multilaminate theory, the numerical integral from a mathematical function is determined by spreading in a sphere area with radius of one. Such a mathematical function can express the changes in the physical properties of a sphere area. For determining the numerical integral, the hypothetical sphere area with a radius of one can be approximated with several flat planes tangential to different points of the sphere area (Figure 2(a)). Each plane has a contact point with the sphere surface; thus, by limiting such planes, the number of contact points or basic points can be defined. When calculating numerical integral, the quantity spread on sphere surface area can be obtained in the mentioned points. The numerical integral is obtained from the continuous function of f(x, y, z) on the sphere surface (the sum of the values of f(x, y, z) in sample points that is multiplied by the weight coefficients related to the points). The number of sample points should be increased in order to reduce errors. The following relation shows the relationship between the numerical integral and the normal integral:

$${\displaystyle {\iint }_{\Omega }f(x,y,z)dxdydz=4\pi {\displaystyle \sum _{i=1}^n{w}_{(i)}f(x_i,y_i,z_i)}}$$ (1)

Ω = sphere area

n = number of points

w_(i) = weight coefficient of point i

_^fi (_^xi , _^yi , _^zi ) = value of function f in point i

Figure 2:
The position of 17 planes a) on the sphere surface b) inside the cubic.

Thus, if the sliding, opening, and closing of each plane are determined by constitutive equations, the sum of the sliding, opening, and closing forms the internal mechanism of the material movement for one point. The total effects of movement or deformation at one point can be obtained by integral addition.

To display the geometry of planes properly, in the present model, the tangential planes for such micro-planes are shown in hemispheres with a radius of one, as shown in Figure 2(a). Figure 2(b) shows their locations inside the unit cubic, while Table 1 indicates the directional cosines, the shear stress directions $${\tau }_{n_1}$$ and $${\tau }_{n_2}$$, and the weight coefficients of 17 planes.

Table1:
Directional cosines and weight coefficients of 17 planes.

Plane(P)	l	m	n	l'	m'	n'	l''	m''	n''	w(P)
1	0.57735	0.57735	0.57735	0.70711	-0.7071	0	0.40825	0.40825	-0.8165	0.0203
2	0.57735	-0.5774	0.57735	0.70711	0.70711	0	0.40825	-0.40825	-0.8165	0.0203
3	-0.5774	0.57735	0.57735	0.70711	0.70711	0	0.40825	-0.40825	0.8165	0.0203
4	-0.5774	-0.5774	0.57735	0.70711	-0.7071	0	0.40825	0.40825	0.8165	0.0203
5	0.70711	0.70711	0	0.70711	-0.7071	0	0	0	1	0.0581
6	-0.7071	0.70711	0	0.70711	0.70711	0	0	0	1	0.0581
7	0.70711	0	0.70711	0	1	0	0.70711	0	-0.70711	0.0301
8	-0.7071	0	0.70711	0	1	0	-0.7071	0	-0.70711	0.0301
9	0	-0.7071	0.70711	1	0	0	0	0.70711	0.70711	0.0301
10	0	0.70711	0.70711	1	0	0	0	-0.70711	0.70711	0.0301
11	1	0	0	0	1	0	0	0	1	0.0383
12	0	1	0	1	0	0	0	0	1	0.0383
13	0	0	1	0	1	0	1	0	0	0.0203
14	0.40825	0.40825	0.8165	0.70711	-0.7071	0	0.57735	0.57735	-0.57735	0.0191
15	0.40825	-0.4082	0.8165	0.70711	0.70711	0	0.57735	-0.57735	-0.57735	0.0191
16	-0.4082	0.40825	0.8165	0.70711	0.70711	0	0.57735	-0.57735	0.57735	0.0191
17	-0.4082	-0.4082	0.8165	0.70711	-0.7071	0	0.57735	0.57735	0.57735	0.0191

If l, m, and n are directional cosines perpendicular to the plane, l’,m’, and n’ are directional cosines for shear stress direction $${\tau }_{n_1}$$, and l”, m”, and n” are directional cosines of shear stress direction $${\tau }_{n_2}$$. According to the results of Sadrnejad and Labibzadeh (2006), stress values in planes can be obtained by the following relation, using matrix algebraic relations and regarding the amount of stress at the point:

$${\text{σ}}_{(P)}=\left[\begin{array}{c}{\sigma }_n\\ {\tau }_{n_1}\\ {\tau }_{n_2}\end{array}\right]=\left[\begin{array}{ccc}l& m& n\\ l^{\prime }& m^{\prime }& n^{\prime }\\ l^{″}& m^{″}& n^{″}\end{array}\right]\left[\begin{array}{ccc}{\sigma }_{xx}& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\sigma }_{yy}& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\sigma }_{zz}\end{array}\right]\left[\begin{array}{c}l\\ m\\ n\end{array}\right]$$ (2)

Here, _^σn is stress perpendicular to the plane, and $${\tau }_{n_1}$$ and $${\tau }_{n_2}$$ are shear stresses in planes in directions 1 and 2 respectively. ^{_σxx , _σyy} , and _^σzz are vertical stresses at the point, while _^τxy , _^τxz , _^τyz ,_^τyx , _^τzx , and _^τzy are shear stresses at the point. In the present research, given the lack of body forces, we have _{^{τxy = τyx, τxz = τzx,}} and _{^{τyz = τzy.}}

After the calculation of stress within planes by the use of transfer relation, elastic and plastic strains are calculated in the planes by the use of relations of the constitutive model. As a result, the strains calculated in planes are transferred to the point by the following relation. The relation is equivalent to the relation of stress transfer from planes to points, as presented by Sadrnejad and Labibzadeh (2006):

$${\epsilon }_{ij}=6\times {\displaystyle \sum _{P=1}^{17}{w}_{(P)}\left[{\text{T}}_{(P)}\right]}\left[{\text{ε}}_{(P)}\right]$$ (3)

ε_(P) is strain that is equivalent to stresses in planes and w _{( P )} is the weight coefficient of planes and ^_εij is the six-component strain in point:

$${\epsilon }_{ij}=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ {\gamma }_{xy}\\ {\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}$$ (4)

T_(P) is plane transfer matrix. The value of T _{( P )} for each plane is calculated from following relation:

$${\text{T}}_{(P)}=\left[\begin{array}{ccc}{l}^2& l{l}^{\prime }& l^{″}l\\ m^2& m^{\prime }m& m^{″}m\\ n^2& n^{\prime }n& n^{″}n\\ lm& 1/2(l^{\prime }m+m^{\prime }l)& 1/2(m^{″}l+l^{″}m)\\ mn& 1/2(m^{\prime }n+n^{\prime }m)& 1/2(m^{″}n+n^{″}m)\\ \ln & 1/2(l^{\prime }n+n^{\prime }l)& 1/2(l^{″}n+n^{″}l)\end{array}\right]$$ (5)

3 CONSTITUTIVE MODEL IN THE PLANE

Among the different models representing sand behaviour in liquefaction soil, the model given by Dafalias and Manzari (2004) has remarkable features. Consistent with the principles of limited soil mechanic, definition of bounding surface, and dilatancy surface, and their dependency on state parameter ψ, and them change during loading and eventually their match with critical surface in failure condition, cause correct prediction of dilation and contraction of dense and loose sand particles (compared to critical state line) and their softening. In the present research, according to the multilaminate framework for the prediction of fabric anisotropy effects, parameters related to the fabric effect of sandy soil particles were omitted from the equation, while low constants were used in different conditions compared to the initial model. According to these features, for the better performance of the mentioned model and to use it in different loading conditions and increase efficiency, the model given by Dafalias and Manzari (2004) was used as the constitutive model in 17 planes via the multilaminate method. In the following section, governing equations to planes and modifications were dealt with using the equations given by Dafalias and Manzari (2004). It should be noted that in all equations, the index ‘P’ in parentheses means the number of the planes and bold expressions show tensor variables.

3.1 Elastic Constitutive Model in the Plane

In the framework of hypo-elastic theory, the elastic shear modulus G _{( P )} and the elastic bulk modulus K _{( P )}, as functions of void ratio e and stress perpendicular to the plane _{^σn ( P )}, are defined as follows:

$$G_{(P)}=G_{0P}{P}_{atm}{(2.97-e)}^2/(1+e){(\frac{{\sigma }_{n(P)}}{P_{atm}})}^{1/2}$$ (6)

and

$$K_{(P)}=\frac{2(1+{\nu }_P)}{3(1-2{\nu }_P)}{G}_{(P)}$$ (7)

Here, _^νP is Poisson's ratio of soil, G _{0 P} is model constant in the plane, and _^Patm is atmospheric pressure.

3.2 Critical State Line

In this model, the soil limit state theory is used in planes. According to this theory, the critical void ratio _{^ec ( P )} and the normal stress _{^σn ( P )} are exponentially inter-related in planes via the following relation:

$$e_{c(P)}=e_{0P}-{\lambda }_{cP}{(\frac{{\sigma }_{n(P)}}{P_{atm}})}^{{\xi }_P}$$ (8)

e_0P and λ_cP are the model constants in the plane. The state parameter ψ _{( P )} is defined as follows for each plane:

$${\psi }_{(P)}=e_{(P)}-e_{c(P)}$$ (9)

e_(P) is equal to e at the beginning of the loading. Bounding, dilatancy, and critical surfaces are defined using this parameter. Depending on the change in this parameter during loading, the aforementioned surfaces were variable as well, leading to the prediction of some features of sand, such as phase change of dilation, softening, and compatibility with limit state theory.

3.3 Yield Surface

Yield surface function in the plane is defined by the following equation:

$$f_{(P)}={[({\text{S}}_{(P)}-{\sigma }_{n(P)}{\text{α}}_{(P)}):({\text{S}}_{(P)}-{\sigma }_{n(P)}{\text{α}}_{(P)})]}^{1/2}-\sqrt{2/3}{\sigma }_{n(P)}{m}_P=0$$ (10)

α_(P) is the deviatoric back stress ratio tensor. This component is used to determine the yield surface axis. The sign ‘:’ means the trace and multiplication of two tensors. Equation f _{( P )} is cone geometry in the deviatoric stress space of the plane. Figure 3 indicates the schematic of yield, dilatancy, bounding, and critical surfaces at the point, and they are observed similarly in the plane.^_mp is one of the plane's constants and the deviatoric stress tensor S _{( P )} in the plane, is calculated as follows:

$${\text{S}}_{(P)}={\text{σ}}_{(P)}-{\sigma }_{n(P)}{\text{I}}_{(P)}$$ (11)

where:

$${\text{I}}_{(P)}=\left[\begin{array}{ccc}\text{1}& \text{0}& \text{0}\\ \text{0}& \text{0}& \text{0}\\ \text{0}& \text{0}& \text{0}\end{array}\right]$$ (12)

r_(P) is stress ratio tensor in the plane. It is calculated by following equation:

$${\text{r}}_{(P)}={\text{S}}_{(P)}/{\sigma }_{n(P)}$$ (13)

Figure 3:
Schematic of yield surface, dilatancy surface, bounding surface, and critical surface at the point and in axes space _^σn ,_{^τn 1} and _{^τn 2} in the plane.

According to equation of f _{( P )}, the yield surface gradient in space r _{( P )} is obtained as follows:

$${\text{L}}_{(P)}=\partial f_{(P)}/\partial {\text{σ}}_{(P)}$$ (14)

n_(P), i.e. the tensor perpendicular to the yield surface in the deviatoric stress space in the plane, is defined as follows:

$${\text{n}}_{(P)}=\frac{{\text{r}}_{(P)}-{\text{α}}_{(P)}}{\Vert {\text{r}}_{(P)}-{\text{α}}_{(P)}\Vert }$$ (15)

where ║║ shows tensor norm.

3.4 Changes in Plastic Strain

Plastic strain increment is calculated from the following relation:

$$d{\text{ε}}_{(P)}^P=〈L_{(P)}〉{\text{R}}_{(P)}$$ (16)

Here, L _{( P )} is the loading index or plastic coefficient, 〈 〉 are Macaulay’s brackets, and 〈x〉 is equal to x if x > 0, while it is equal to zero if x ≤ 0. R _{( P )} is the vector direction $$d{\text{ε}}_{(P)}^P$$. In this model and in general, R _{( P )} ≠ ∂f _{( P )} / ∂ σ _{( P )} meaning that the non-associated flow rule is governed in plasticity R _{( P )} is divided into two deviatoric and volumetric parts using the following relation. The equation is as follows:

$${\text{R}}_{(P)}={\text{R}}_{(P)}^{\prime }+\frac{1}{3}{D}_{(P)}{\text{I}}_{(P)}=B_{(P)}{\text{n}}_{(P)}+C_{(P)}({\text{n}}_{(P)}^2-\frac{1}{3}{\text{I}}_{(P)})+\frac{1}{3}{D}_{(P)}{\text{I}}_{(P)}$$ (17)

Here, R’ _{( P )} is the deviatoric part R _{( P )} B _{( P )} and C _{( P )} are used to study the effect of Lode's angle on the plane, as follows:

$$B_{(P)}=1+\frac{3}{2}\times \frac{(1-C_P)}{C_P}\times g_{(P)}\times \cos 3{\theta }_{(P)}$$ (18)

$$C_{(P)}=3\sqrt{\frac{3}{2}}\times \frac{(1-C_P)}{C_P}\times g_{(P)}$$ (19)

C_(P) is the constant in the plane and θ _{( P )} is the effect of Lode's angle in the plane. It is defined as follows:

$$\cos 3{\theta }_{(P)}=\sqrt{6}trc{\text{n}}_{(P)}^3$$ (20)

g_(P) is the interpolation function to determine bounding, dilatancy, and critical surfaces in the different conditions of stress path. It is defined by following relation:

$$g_{(P)}=\frac{2C_P}{(1+C_P)-(1-C_P)COS3{\theta }_{(P)}}$$ (21)

The dilation parameter for the plane D _{( P )} is defined by the following relation:

$$D_{(P)}=A_{d(P)}({\text{α}}_{\theta (P)}^d-{\text{α}}_{(P)}):{\text{n}}_{(P)}$$ (22)

The dilatancy surface $${\text{α}}_{\theta (P)}^d$$ is defined as follows:

$${\text{α}}_{\theta (P)}^d=\sqrt{2/3}\times [g_{(P)}\times M_P\times \exp (n_P^d\times {\psi }_{(P)})-m_P]{\text{n}}_{(P)}$$ (23)

In this relation, $$n_P^d$$ is the model constant in the plane. Dafalias and Manzari (2004) considered _^Ad parameter (equivalent to _{^Ad ( P )} in the plane) to predict the fabric effect. The value of _{^Ad ( P )} is calculated by the following relation:

$$A_{d(P)}=A_{0P}(1+〈{\text{Z}}_{(P)}:{\text{n}}_{(P)}〉)$$ (24)

A_0P is the model constant in the plane. The fabric dilatancy tensor rate in the plane is calculated as follows:

$$d{\text{Z}}_{(P)}=-C_{ZP}〈-d{\epsilon }_{{}_{V(P)}}^P〉(Z_{\max P}{\text{n}}_{(P)}+{\text{Z}}_{(P)})$$ (25)

Given the elimination of _^CZP and _{^{Zmax P}} parameters of fabric tensor, or assuming a zero value for them, in the multilaminate proposed method, modified equation D _{( P )} in multilaminate method is as follows:

$$D_{(P)}=A_{0P}({\text{α}}_{\theta (P)}^d-{\text{α}}_{(P)}):{\text{n}}_{(P)}$$ (26)

The increment of plastic volumetric strain in the plane is calculated as follows:

$$d{\epsilon }_{{}_{V(P)}}^P=〈L_{(P)}〉D_{(P)}$$ (27)

The rate of back stress ratio tensor dα _{( P )} is calculated as follows:

$$d{\text{α}}_{(P)}=〈L_{(P)}〉2/3h_{(P)}({\text{α}}_{\theta (P)}^b-{\text{α}}_{(P)})$$ (28)

$${\text{α}}_{\theta (P)}^b$$ and h _{( P )} are calculated by the following relations:

$${\text{α}}_{\theta (P)}^b=\sqrt{2/3}\times [g_{(P)}\times M_P\times \exp (-n_P^b\times {\psi }_{(P)})-m_P]{\text{n}}_{(P)}$$ (29)

$$h_{(P)}=\frac{b_{0(P)}}{({\text{α}}_{(P)}-{\text{α}}_{ini(P)}):{\text{n}}_{(P)}}$$ (30)

_^αini(P) is the value ofα _{( P )} at the start of each new step of loading and unloading. Coefficient b _{0( P )} is calculated using the following equation:

$$b_{0(P)}=G_{0P}{h}_{0P}(1-C_{hP}e){(\frac{{\sigma }_{n(P)}}{P_{atm}})}^{-1/2}$$ (31)

In this equation, _^ChP and h _{0 P} are the model constants in the planes. According to mathematical calculations, the loading index L _{( P )} is estimated by the following relation:

$$\begin{array}{l}{L}_{(P)}=(2\times G_{(P)}{\text{n}}_{(P)}:d{\text{e}}_{(P)}-{\text{n}}_{(P):}:{\text{r}}_{(P)}d{\epsilon }_{V(P)}\times K_{(P)})/(K_{P(P)}+\\ 2G_{(P)}(B_{(P)}+C_{(P)}tr{\text{n}}_{(P)}^3)-K_{(P)}{D}_{(P)}{\text{n}}_{(P)}:{\text{r}}_{(P)})\end{array}$$ (32)

The plastic coefficient _{^KP ( P )} is calculated as follows:

$$K_{P(P)}=2/3\times {\sigma }_{n(P)}\times h_{(P)}({\text{α}}_{\theta (P)}^b-{\text{α}}_{(P)}):{\text{n}}_{(P)}$$ (33)

Finally, given that $$d{\epsilon }_{V(P)}=d{\epsilon }_{V(P)}^e+d{\epsilon }_{V(P)}^P$$ and $$d{\text{e}}_{(P)}=d{\text{e}}_{(P)}^e+d{\text{e}}_{(P)}^P$$, the stress increment in the plane dσ _{( P )} is calculated by the following relation:

$$\begin{array}{l}d{\text{σ}}_{(P)}=2G_{(P)}d{\text{e}}_{(P)}^e+K_{(P)}d{\epsilon }_{V(P)}{\text{I}}_{(P)}-〈L_{(P)}〉\{2G_{(P)}[B_{(P)}{\text{n}}_{(P)}+\\ C_{(P)}({\text{n}}_{(P)}^2-1/3{\text{I}}_{(P)})]+K_{(P)}{D}_{(P)}{\text{I}}_{(P)}\}\end{array}$$ (34)

where, the elastic deviatoric strain $$d{\text{e}}_{(P)}^e$$ is calculated according to the following equation:

$$d{\text{e}}_{(P)}^e=d{\text{s}}_{(P)}/(2*G_{(P)})$$ (35)

3.5 Undrained Condition

In the present research, excess pore water pressure is calculated to model the undrained condition using Naylor method (Naylor, 1974; Tam, 1992):

$$\delta \text{u}=K_u{\text{m}}^T\delta {\epsilon }_v$$ (36)

Here, _^δεv is the increment of the total volumetric strain and volumetric modulus of soil and water_^Ku , and _^mT are calculated by the following relations:

$$\frac{1}{K_u}=\frac{(1-n)}{K_S}+\frac{n}{K_W}$$ (37)

$${\text{m}}^T=\left[\begin{array}{cccccc}1& 1& 1& 0& 0& 0\end{array}\right]$$ (38)

In these relations, n is the soil porosity, _^KW is the water volumetric modulus, and _^KS is the volumetric modulus of soil granules.

4 DETERMINING MODEL CONSTANTS

The model constants given by Dafalias and Manzari (2004) in previous research works have been obtained for Toyoura sand, based on the experiment results (Taiebat and Dafalias, 2008; Taiebat et al., 2010). In the present research, given similar constants in all planes, sensitivity analysis is done for the undrained condition with p ₀ = 1000kPa and e ₀ = 0.833 as mean values of pressure and void ratio; their values are shown in Table 2.

Table2:
Values of constants in the plane.

Constant	Parameter	Value
Elasticity	G 0P	120
	ν P	0.35
Critical State	M P	0.72
	C P	0.93
	λCP	0.019
	e 0 P	0.955
	ξP	0.85
Yield Surface	m P	0.02
Plastic Modulus	h oP	4
	C hP	1
	$$n_P^b$$	1.1
Dilatancy	A 0 P	0.5
	$$n_P^d$$	4

5 RESULTS AND DISCUSSION

In the present research, the fabric effect has been modelled using the multilaminate model framework. Based on microscopic observations, the sandy soil particles should demonstrate contraction behaviour during inverse load, when they are placed in dilatancy condition. To modify their model, Manzari and Dafalias (1997) used the fabric effect in the equations of factor D, using constants _^CZ and Z _max Therefore, sand modelling in cyclic loading was correctly predicted and effective zero stress and soil liquefaction were calculated. Figure 4 shows the results of cyclic loading without the fabric effect parameters (Manzari and Dafalias, 1997). Figure 5 indicates cyclic loading with the fabric effect parameters (Dafalias and Manzari, 2004). According to the ability of the model given by Dafalias and Manzari (2004) in different monotonic and cyclic loadings and modelling of drained and undrained conditions, the fabric effect was considered in this research by another method, which is able to predict such behaviour with lower constants. Therefore, by removing the constants related to the fabric effect, numerical modelling was done in the multilaminate framework. Figure 6 shows the results of this model, with 13 constants compared to the 15 soil constants in Dafalias and Manzari's (2004) model. It is observed that multilaminate numerical framework is able to model the fabric effect in cyclic condition. Therefore, a more simplified model is presented, due to its remarkable features such as different loading conditions, prediction of softening, and limit state condition.

Figure 4:
The results of cyclic loading (Dafalias-Manzari method) without the effect of fabric constants at e ₀ = 0.808 and p ₀ = 294kPa.

Figure 5:
The results of cyclic loading (Dafalias-Manzari method) with effect of fabric constants at e ₀ = 0.808 and p ₀ = 294kPa.

Figure 6:
The results of cyclic loading via multilaminate method and removal of fabric constants.

To study more about the effect of the removal of constants related to the fabric in the multilaminate framework in different loading conditions, void ratio, and confining pressure, the numerical results of constants _^CZP = 600 and Z _{max P} = 4 (resulting from calibration) were compared to their removal conditions or zero constants (Figures 7-12). To verify the numerical simulation of the model, the results were compared to the results obtained from the triaxial standard tests of Verdugo and Ishihara (1996). According to the wide changes of confining pressure on the samples from 100kPa to 3000kPa and void ratio from 0.735 to 0.907, results were evaluated in various void and confining pressure conditions during loading and unloading.

Figure 7:
Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in drained standard triaxial test at p ₀ = 100kPa with and without fabric constants.

Figure 8:
Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in drained standard triaxial test at p ₀ = 500kPa with and without fabric constants.

Figure 9:
Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in the undrained standard triaxial test at e ₀ = 0.735 with and without fabric constants.

Figure 10:
Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in the undrained standard triaxial test at e ₀ = 0.833 with and without fabric constants.

Figure 11:
Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in the undrained standard triaxial test at e ₀ = 0.907 with and without fabric constants.

Figure 12:
The comparison of experimental (data from Dafalias and Manzari (2004)) results (a, b) and multilaminate modelling for cyclic loading standard triaxial test at e ₀ = 0.808 and p ₀ = 294kPa, with and without fabric constants(c,d).

The drained condition, with p ₀ = 100kPa and different e ₀ for the multilaminate model, has been compared to the experimental results in Figure 7. This comparison was done with the above condition for p ₀ = 500 kPa, as shown in Figure 8. The results show the prediction ability of the model in different conditions, such as softening of dense sand particles.

Figure 9 indicate the results for undrained condition, with e ₀ = 0.735 in different confining pressures. An important point in the figures is that the steeper gradient of curve p-q is for low confining pressures (p ₀ = 100kPa), which matches with the theoretical discussions suggesting the increase in soil angle of internal friction at low pressures and its reduction at high confining pressures.

In Figure 10, calibration of the results of the plane constants was done with this void ratio (e ₀ = 0.833); thus, it is seen to be more in agreement with the experimental results.

The ability of this model to predict liquefaction appropriately in loose sandy soil with void ratio e ₀ = 0.907 is shown in Figure 11.

Finally, Figure 12 shows the ability of the model in cyclic loading. Among the remarkable features of the model is the prediction of inverse loading in planes and its effect on the point. The reason for different results in cyclic condition is that cyclic test was done 20 years ago, compared to monotonic tests. This issue can be related by testing method and instrumentation (Dafalias and Manzari, 2004).

As seen in Figures 7-12, there is not much difference after the removal of parameters related to the fabric effect using the multilaminate method. This can be due to fundamental changes in the calculational principles of the multilaminate method and material behavioural features in the point, resulting from 17 planes instead of three planes. The closeness between the results of the multilaminate model and the experimental results proves the model capability.

The induced anisotropy effect on post-liquefaction sand particles in 17 directions was studied using related Figures (13-16) in the standard triaxial test in monotonic and cyclic conditions. Compression and extension were assumed positive and negative respectively in all the modelling. To this end, a sand particle with p ₀ = 100kPa and e ₀ = 907 in monotonic and undrained loading and unloading conditions was used for the modelling.

Figure 13:
The shear stress versus normal stress on 17 planes in monotonic undrained condition, with p ₀ = 1000kPa and e ₀ = 0.907.

Figure 14:
The shear stress versus shear strain on 17 planes in monotonic undrained condition, with p ₀ = 1000kPa and e ₀ = 0.907.

Figure 15:
The shear stress versus normal stress on 17 planes in cyclic undrained condition, with p ₀ = 294kPa and e ₀ = 0.808.

Figure 16:
The shear stress versus shear strain on 17 planes in cyclic undrained condition, with p ₀ = 294kPa and e ₀ = 0.808.

Figure 13 shows the stress path on 17 planes. As seen in Planes 9, 10, 11, 12, and 13, no shear stress is present; there is only pressure. Therefore, such planes are inactive shear and there is no shear strain in these planes (Figure 14). Therefore, these planes have no contribution in the post-liquefaction lateral spreading. Planes 5, 6, 7, and 8, Planes 1, 2, 3, and 4, and Planes 14, 15, 16, and 17 all follow a stress path to reach failure and liquefaction. This can be studied in Figure 14, where Planes 3 and 4 experience more strain compared to Planes 1 and 2 due to the effect of the similar stress path. The resultant strain determines the movement of soil particles and displacement created in the liquefaction soil in different directions after liquefaction.

Eventually, Figures 15 and 16 evaluated the soil behaviour in 17 planes in undrained and cyclic condition under the loading effect of the triaxial standard test with p ₀ = 294kPa and e ₀ = 0.808.

Figure 15 shows the stress path in the planes. As observed, the value of cyclic shear stress is maximum in Planes 5 and 6 and minimum in Planes 14, 15, 16, and 17. The stress path on each plane follows a different pattern to reach liquefaction. It is noteworthy that Planes 7, 8, 9, 10, 11, 12, and 13 are inactive in cyclic conditions in the triaxial standard test, and no shear stress is produced.

Figure 16 indicates shear strain versus shear stress in 17 planes. It should be noted that more shear strain is produced in some planes, such as Planes 5 and 6, while Planes 16 and 17 experience less strain. Therefore, the direction of the post-liquefaction strain in cyclic condition is determined for planes that have more contribution to strain production.

Figures 13-16 show that the status of planes, their strain history, and the values of shear strength and strain produced in them in different loading conditions confirm the ability of the model in predicting the induced anisotropy and post-liquefaction occurrences.

6 CONCLUSION

In the present research, in addition to the study effect of constant reductions in an advanced constitutive model in the multilaminate framework, the effect of induced anisotropy in pre- and post-liquefaction sandy soil is evaluated. In this framework, the constitutive model in the planes is simplified and the costs of experiments and calculations are reduced in practice. In the multilaminate method, the constitutive model was used in 17 planes with different directions instead of in three planes. This can increase the accuracy of the calculations. If the conditions are provided to define different constants in various planes via the experimental results based on soil properties under the initial condition, or even if different constitutive equations are used for different planes, the effect of initial anisotropic condition on soil behaviour can be properly evaluated after loading with this method.

The ability of the multilaminate model in predicting sand particle behaviour in different conditions and the start of liquefaction for different values of pressure and void ratio, boundary condition of drained and undrained, cyclic, and monotonic loadings was compared to the experimental results of triaxial standard test and good results were obtained.

Stress-strain behaviour was studied in 17 planes to investigate the model performance in induced anisotropy, liquefaction occurrence, and post-liquefaction conditions. Different values of stress and strain in planes show that the induced anisotropy is affected by loading condition and plane arrangement; it proves the ability of the new model for the estimation of such anisotropy. This is especially important in liquefaction soil under different topographic condition, in situ stress and strain, to predict their behaviour and displacements occurring after liquefaction.

Multilaminate model framework provides rotation of the principal axis of stress and lack of co-axiality of principal stress and strain and involves the effects of partial changes in boundary condition, as it protects directional effects in material behaviour. Such subjects and other capabilities of the model can be evaluated in further research works.

References

• Basari, E., Ozden, G. (2013). Post-liquefaction volume change in micaceous sandy soils of Old Gediz River Delta , Acta Geotechnica Slovenica 1:33-40.
• Bazant, Z. P., Xiang, Y., Adley, M. D., Prat, P. C., Akers, S. A. (1996). Microplane model for concrete: II: data delocalization and verification, Journal of Engineering Mechanics;122(3):255-262.
• Caner, F. C., Bazant, Z. P. (2000). Microplane model M4 for concrete. II: Algorithm and calibration, Journal of Engineering Mechanics 126(9):954-961.
• Chang, C. S., Hicher, P. Y. (2005). An elasto-plastic model for granular materials with microstructural consideration, International Journal of Solids and Structures 42:4258-4277.
• Dafalias, Y. F. (1986). Bounding surface plasticity. I: Mathematical foundation and hypoplasticity, Journal of Engineering Mechanics 112(9):966-987.
• Dafalias, Y. F., Manzari, M. T. (2004). Simple Plasticity Sand Model Accounting for Fabric Change Effects, Journal of Engineering Mechanics 130:622-634.
• Dafalias, Y. F., Papadimitriou, A. G., Li, X. S. (2004). Sand Plasticity Model Accounting for Inherent Fabric Anisotropy, Journal of Engineering Mechanics 130(11):1319-1333.
• Elgmal, A., Yang, Z., Parra, E. (2002). Computational modeling of cyclic mobility and post-liquefaction site response, Soil Dynamics and Earthquake Engineering 22:259-271.
• Gao, Z., Zhao, J., Li, X. S., Dafalias, Y. F. (2014). A critical state sand plasticity model accounting for fabric evolution, International Journal for Numerical and Analytical Methods in Geomechanics 38(4):370-390.
• Ghadrdan, M., Sadrnejad, S. A. (2015). Numerical evaluation of geomaterials behavior upon multiplane damage Model, Computers and Geotechnics 68:1-7.
• Guo, N., Zhao, J. (2013). The signature of shear-induced anisotropy in granular media, Computers and Geotechnics 47:1-15.
• Guo, P. (2014). Coupled effects of capillary suction and fabric on the strength of moist granular materials, Acta Mechanica 225:2261-2275.
• Hareb, H., Doanh, T. (2012). Probing into the strain induced anisotropy of Hostun RF loose sand, Granular Matter 14:589-605.
• Hicher, P. Y., Chang, C. S. (2007). A microstructural elastoplastic model for unsaturated granular materials, International Journal of Solids and Structures 44:2304-2323.
• Ishikawa, A., Zhou, Y. G., Shamoto, Y., Mano, H., Chen, Y. M., Ling, D. S. (2015). Observation of post- liquefaction progressive failure of shallow foundation in centrifuge model tests, Soils and Foundations 55(6):1501-1511.
• Khalili, Y., Mahboubi, A. (2014). Discrete simulation and micromechanical analysis of two-dimensional saturated granular media, Particuology 15:38-150.
• Kruyt, N. P., Rothenburg, L. (2016) A micromechanical study of dilatancy of granular materials, Journal of the Mechanics and Physics of Solids 95:411-427.
• Kruyt, N. P., Rothenburg, L. (2014). On micromechanical characteristics of the critical state of two-dimensional granular materials, Acta Mechanica 225:2301-2318.
• Lashkari, A. (2009). A constitutive model for sand liquefaction under rotational shear, Iranian Journal of Science and Technology 33:31-48.
• Li, X. S. (2002). A sand model with state dependent dilatancy, Geotechnique 52(3):173-186.
• Li, X., Yu, H. S. (2014). Fabric, force and strength anisotropies in granular materials: a micromechanical insight, Acta Mechanica 225:2345-2362.
• Liu, Y., Chang, C. S., Wu, S. C. (2016). A Simple One-Scale Constitutive Model for Static Liquefaction of Sand-Silt Mixtures, Latin American Journal of Solids and Structures 13:2190-2218.
• Louis, L., Baud, P., Wong, T. F. (2009). Microstructural Inhomogeneity and Mechanical Anisotropy Associated with Bedding in Rothbach Sandstone, pure and applied geophysics 66:1063-1087.
• Manzari, M. T., Dafalias, Y. F. (1997). A critical state two-surface plasticity model for sands, Geotechnique 47(2):255-272.
• Naylor, D. J. (1974). Stresses in nearly incompressible materials by finite elements with application to the calculation of excess pore pressures, International Journal for Numerical Methods in Engineering 8(3):443-460.
• Sadrnejad, S. A. (2007). A general multi-plane model for post-liquefaction of sand, Iranian Journal of Science & Technology, Transaction B, Engineering 31(B1):123-141.
• Sadrnejad, S. A. (2011). Soil Plasticity and Modeling (2th ed). Khajeh Nasir Toosi University of Technology(Tehran).
• Sadrnejad, S. A. (2009). Semi-Micro Bounding Surface Model for Anisotropic Sand, Electronic Journal of Geotechnical Engineering 14:1-25.
• Sadrnejad, S. A., Daryan, A. S., Ziaei, M. (2009). A Constitutive Model for Multi-Line Simulation of Granular Material Behavior Using Multi-Plane Pattern, Journal of Computer Science 5(11):822-830.
• Sadrnejad, S. A., Karimpour, H. (2011). Drained and undrained sand behaviour by multilaminate bounding surface model, International Journal of Civil Engineering 9(2):111-125.
• Sadrnejad, S. A., Labibzadeh, M. (2006). A Continuum/Discontinuum Micro Plane Damage Model for Concrete. International Journal of Civil Engineering 4(4):296-313.
• Shamoto, Y., Zhang, J. M., Tokimatsu, K. (1998). New charts for predicting large residual post-liquefaction ground deformation, Soil Dynamics and Earthquake Engineering 17 :427-438.
• Taiebat, M., Jeremic, B., Dafalias, Y. F., Kaynia, A. M., Cheng, Z. (2010). Propagation of seismic waves through liquefied soils, Soil Dynamics and Earthquake Engineering 30:236-257.
• Taiebat, M., Dafalias, Y. F. (2008). SANISAND: Simple anisotropic sand plasticity model, International Journal for Numerical and Analytical Methods in Geomechanics 32:915-948.
• Tam, H. K. (1992). Some applications of cam-clay in numerical analysis. Ph.D. Thesis, City University London, England.
• Tang Tron Tran, H., Wong, H., Dubujet, Ph, Doanh, T. (2014). Simulating the effects of induced anisotropy on liquefaction potential using a new constitutive model, International Journal for Numerical and Analytical Methods in Geomechanics 38(10):1013-1035.
• Verdugo, R., Ishihara, K. (1996). The steady state of sandy soils, Soils and Foundations 36(2):81-91.
• Voyiadjis, G. Z., Thiagarajan, G., Petrakis, E. (1995). Constitutive modelling for granular media using an anisotropic distortional yield model, Acta Mechanica 110:151-171.
• Wang, G., Xie, Y. (2014). Modified Bounding Surface Hypoplasticity Model for Sands under Cyclic Loading, Journal of Engineering Mechanics 140:91-101.
• Wang, R., Zang, J. M., Wang, G. (2014). A unified plasticity model for large post-liquefaction shear deformation of sand, Computers and Geotechnics 59:54-66.
• Wang, S., Luna, R., Onyejekwe, S. (2015). Postliquefaction behavior of low-plasticity silt at various degrees of reconsolidation, Soil Dynamics and Earthquake Engineering 75:259-264.
• Wang, Z. L., Dafalias, Y. F., Shen, C. K. (1990). Bounding surface hypoplasticity model for sand, Journal of Engineering Mechanics 116(5):983-1001.
• Yan, W. M. (2009). Fabric evolution in a numerical direct shear test, Computers and Geotechnics 36:597-603.
• Yang, Z. X., Li, X. S., Yang, J. (2008). Quantifying and modeling fabric anisotropy of granular soils, Geotechnique 58(4):237-248.
• Yang, Z. X., Yang, J., Wang, L. Z. (2013) Micro-scale modeling of anisotropy effects on undrained behavior of granular soils, Granular Matter 15:557-572.
• Ye, B., Ye, G., Zang, F. (2012). Numerical modeling of changes in anisotropy during liquefaction using a generalized constitutive model, Computers and Geotechnics 42:62-72.
• Yu, H., Zeng, X., Li, B., Ming, H. (2013). Effect fabric anisotropy on liquefaction of sand, Journal of Geotechnical and Geoenvironmental Engineering 139:765-774.
• Zhang, J. M., Wang, G. (2012). Large post-liquefaction deformation of sand, part I: physical mechanism, constitutive description and numerical algorithm, Acta Geotechnica 7:69-113.
• Zhao, J., Guo, N. (2013). Unique critical state characteristics in granular media considering fabric anisotropy, Geotechnique 63(8):695-704.

Publication Dates

• Publication in this collection
  16 Aug 2017

History

• Received
  19 Mar 2017
• Reviewed
  05 June 2017
• Accepted
  14 June 2017