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

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 preand post-liquefaction behaviour in granular soil. In the framework of the multilaminate method, 17 planes with predetermined 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 preand 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.

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-Latin American Journal of Solids and Structures 14 (2017) 1569-1593 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 postliquefaction 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 postliquefaction 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 Latin American Journal of Solids and Structures 14 (2017) 1569-1593 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 modelalso 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.

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 planesdue 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: 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.
Latin American Journal of Solids and Structures 14 (2017) 1569-1593   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 : 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): T is plane transfer matrix.The value of ( )

P
T for each plane is calculated from following relation: Latin American Journal of Solids and Structures 14 (2017) 1569-1593

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 y , 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.

Elastic Constitutive Model in the Plane
In the framework of hypo-elastic theory, the elastic shear modulus e 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.

Yield Surface
Yield surface function in the plane is defined by the following equation: where: r is stress ratio tensor in the plane.It is calculated by following equation: According to equation of ( ) P f , the yield surface gradient in space ( ) P r is obtained as follows: ( ) ( ) ( ) n , 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.

Changes in Plastic Strain
Plastic strain increment is calculated from the following relation: Here, ( ) P L is the loading index or plastic coefficient, are Macaulay's brackets, and x is equal to x if 0 x > , while it is equal to zero if is the vector direction C are used to study the effect of Lode's angle on the plane, as follows: (1 ) 3 1 c o s 3 2 (1 ) 3 3 2 P C is the constant in the plane and ( ) P q is the effect of Lode's angle in the plane.It is defined as follows: cos 3 6 g is the interpolation function to determine bounding, dilatancy, and critical surfaces in the different conditions of stress path.It is defined by following relation: 1) (1 ) 3 The dilation parameter for the plane ( ) P D is defined by the following relation: ): The dilatancy surface ( ) α is defined as follows: A is calculated by the following relation: (1 : ) 0P A is the model constant in the plane.The fabric dilatancy tensor rate in the plane is calculated as follows: ( ) ( ) max ( ) ( ) ( ) ): The increment of plastic volumetric strain in the plane is calculated as follows: ( ) ( ) ( ) The rate of back stress ratio tensor ( ) P dα is calculated as follows: 2 / 3 ( ) L is estimated by the following relation: dσ is calculated by the following relation:

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 de is the increment of the total volumetric strain and volumetric modulus of soil and water u K , and T m are calculated by the following relations: In these relations, n is the soil porosity, W K is the water volumetric modulus, and S K is the volumetric modulus of soil granules.

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)

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 Z C and max Z .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 600 ZP C = and max 4 P Z = (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.sures.An important point in the figures is that the steeper gradient of curve p-q is for low confining pressures ( 0 100 p k P a =

Latin
), 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 ( 0 0.833 e = ); 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 0 0.907 e = 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 0 1000 p k P a = and 0 0.907 e = 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.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][14][15][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.

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 postliquefaction 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 coaxiality 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.
Modification of a Constitutive Model in the Framework of a Multilaminate Method for Post-Liquefaction Sand 1 INTRODUCTION

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

Figure 2 :
Figure 2: The position of 17 planes a) on the sphere surface b) inside the cubic.
to stresses in planes and ( )

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

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

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

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

e
for the multilaminate model, has been compared to the experimental results in Figure 7.This comparison was done with the above condition for 0 500 p k P a = , 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 7 :Figure 8 :
Figure 7: Comparison of experimental results and multilaminate modelling for monotonic loading and unloading in drained standard triaxial test at 0 100 p k P a = with and without fabric constants.

Figure 9
Figure 9 indicate the results for undrained condition, with 0 0.735 e = in different confining pres-

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

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

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

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

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

Figure 15 :Figure 16 :
Figure 15: The shear stress versus normal stress on 17 planes in cyclic undrained condition, with 0 294 p k P a = and 0 0.808 e = .

Table 1 :
Directional cosines and weight coefficients of 17 planes.
n s is stress perpendicular to the plane, and 1 n t and 2 n t are shear stresses in planes in directions 1 and 2 respectively.xx s , yy s , and zz s are vertical stresses at the point, while xy t , xz t , Latin American Journal of Solids and Structures 14 (2017) 1569-1593 yz t , yx t , zx t , and zy t are shear stresses at the point.In the present research, given the lack of body forces, we have xy In this model, the soil limit state theory is used in planes.According to this theory, the critical void ratio K , as functions of void ratio e and stress perpendicular to the plane ( ) P n is Poisson's ratio of soil, 0P G is model constant in the plane, and atm P is atmospheric pressure.Latin American Journal of Solids and Structures 14 (2017) 1569-1593 3.2 Critical State Line

Table 2 :
Values of constants in the plane.
. In the present research, given similar constants in all planes, sensitivity analysis is done void ratio; their values are shown in Table2.Latin American Journal ofSolids and Structures 14 (2017) 1569-1593