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.

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:

Ω = sphere area

*n* = number of points

*w*
_{(}
_{
i
}
_{)} = weight coefficient of point *i*

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

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

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
*l”*, *m”*, and *n”* are directional cosines of shear stress direction
^{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:

Here, _{
σn
} is stress perpendicular to the plane, and
_{
σ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)}:

**
ε
**

_{(}

_{ 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:

**
T
**

_{(}

_{ P }

_{)}is plane transfer matrix. The value of

*T*_{(}

_{ P }

_{)}for each plane is calculated from following relation:

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:

and

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*
_{0}
_{
P
} and λ_{
cP
} are the model constants in the plane. The state parameter *ψ*
_{(}
_{
P
}
_{)} is defined as follows for each plane:

*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:

**
α
**

_{(}

_{ 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:

where:

**
r
**

_{(}

_{ P }

_{)}is stress ratio tensor in the plane. It is calculated by following equation:

According to equation of *f*
_{(}
_{
P
}
_{)}, the yield surface gradient in space **
r
**

_{(}

_{ P }

_{)}is obtained as follows:

**
n
**

_{(}

_{ P }

_{)}, i.e. the tensor perpendicular to the yield surface in the deviatoric stress space in the plane, is defined as follows:

where ║║ shows tensor norm.

3.4 Changes in Plastic Strain

Plastic strain increment is calculated from the following relation:

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

*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:

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:

*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:

*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:

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

The dilatancy surface

In this relation,
^{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*
_{0}
_{
P
} is the model constant in the plane. The fabric dilatancy tensor rate in the plane is calculated as follows:

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:

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

The rate of back stress ratio tensor **
dα
**

_{(}

_{ P }

_{)}is calculated as follows:

*h*
_{(}
_{
P
}
_{)} are calculated by the following relations:

_{
α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:

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:

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

Finally, given that
**
dσ
**

_{(}

_{ P }

_{)}is calculated by the following relation:

where, the elastic deviatoric strain

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}):

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:

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*
_{0} = 1000*kPa* and *e*
_{0} = 0.833 as mean values of pressure and void ratio; their values are shown in Table 2.

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.

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 100*kPa* to 3000*kPa* and void ratio from 0.735 to 0.907, results were evaluated in various void and confining pressure conditions during loading and unloading.

The drained condition, with *p*
_{0} = 100*kPa* and different *e*
_{0} for the multilaminate model, has been compared to the experimental results in Figure 7. This comparison was done with the above condition for *p*
_{0} = 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} = 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*
_{0} = 100*kPa*), 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} = 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} = 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*
_{0} = 100*kPa* and *e*
_{0} = 907 in monotonic and undrained loading and unloading conditions was used for the modelling.

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*
_{0} = 294*kPa* and *e*
_{0} = 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.