Fabric Behavior of Sand in Post-liquefaction

An anisotropic plasticity model for post-liquefacti on of the undrained behavior of sand is presented. The model incorporates the critical/stea dy state concept that postulates the existence of a state where sand continuously deforms at a certain constant effective stress depending two main parameters of both initial bulk parameters (void ra tio or relative density) and the stress level (mean stress). The local instability of saturated sand wi thin post-liquefaction is highly dependent on the residual inherent/induced anisotropy, bedding plane eff cts and stress/strain path. Most of the models developed using stress/strain invariants are not ca pable of identifying the parameters depending on orientation such as fabric. This is mainly because str ss/strain invariants are quantities similar to scalar quantities and not capable of carrying directional i formation with them. The constitutive equations o f the model are derived within the context of non-lin ear elastic behavior of the whole medium and plasti c sliding of interfaces of predefined multi-planes. T he proposed multi-plane based model is capable of predicting the behavior of soils on the basis of pl astic sliding mechanisms, elastic behavior of parti cles and possibilities to see the micro-fabric effects a s n tural anisotropy as well as induced anisotropy in plasticity. The model is capable of predicting the behavior of soil under different orientation of bed ding plane, history of strain progression during the app lication of any stress/strain paths. The influences of rotation of the direction of principal stress and s train axes and induced anisotropy are included in a rational way without any additional hypotheses. The spatial strength distribution at a location as an approximation of probable mobilized sliding mechani sm is proposed as an ellipsoid function built up on bedding plane.


INTRODUCTION
Saturated deposits composed mainly of cohesion less material have shown to be very likely to undergo flow-type of failure associated with large deformation named post-liquefaction. This type of failure can occur during fast static loading or earthquake. A typical example is during earthquake loading where the strongest perturbation occurs during a few seconds. The tendency of soil mass is to contract due to variation in total confining stress and due to shear induced, resulting instead in an increase of pore water pressure and the reduction of effective stresses. In the extreme case, the effective stresses may become so small (or even zero) resulting in complete loss of strength and the soil may flow in a manner resembling a liquid.
Assuming the validity of continuum mechanics principles, two basic approaches may be considered as one-phase modelling and two-phase modelling. Large deformation features observed in post-liquefaction behavior of sand is simply employed in the first approach, while the second approach is not easily extendable to large deformation. To model the liquefied soil as a one-phase material, the governing rules can be still subdivided into two types as solid mechanic approach and fluid mechanics approach. According to the first approach, the liquefied soil is assumed to be either elastic or elastic-plastic solid having properties of residual state. However, in the second approach, the liquefied soil is treated as a Newtonian fluid having only viscous resistance against shearing and may be modelled using either Navier-Stocks or Euler formulation. However, experimental researches have revealed that the mechanical behavior of liquefied soil is neither solid nor fluid. It seems that it has the characteristics of both solid and fluid like behaviors. These characteristics lead the solid behavioral approach to assume a non-linear elasticity and a calibrated form of plasticity associated with a special flexibility in plasticity modulus.
The steady state concept has been defined as the state under which the mass of cohesionless soil continuously deforms under constant volume, effective stresses and velocity, also varied with void ratio and confining stress. However, despite differences in formal definitions, Pooroshasb [1] argued that the critical state and the steady state concepts are essentially the same and proposed that an appropriate terminology would be the ultimate state. Experimental evidence presented by Been et al. [2] show the critics and steady states are indeed the same. The steady state concept has been used mainly in the study of the instability of sand under large deformation. In order to reliably compute the safety and stability of a cohesionless soil mass against the flow-type failure, it is important to be able to predict the behavior of the soil at all stages of loading, not just at the ultimate state.
Li and Dafalias [3] pointed out that the classical stress ditatancy approach presented by Rowe (1962) in its exact form ignored the extra energy loss due to the static and dynamic constraints at particle contacts and led to a unique relationship between the stress ratio and dilatancy. Furthermore, it has been shown by Manzari and Dafalias [4] , Wan et al. [5] , Li and Dafalias [3] that in order to model sand behavior over a full range of density states, additional dependence of dilatancy on the material internal states and also sliding mechanism are needed and the material state must be described in reference to the critical/steady state line. Accordingly, any local sliding instability may take its limited share in global behavior and through properly controlled scheme there will be the possibility of proceeding deformability of post-liquefied zone.
In addition, due mainly to the process of deposition under earth gravity, the behavior of in situ sand is inherently anisotropy, meaning the stress-strainstrength relations for the same sand may vary as the stress tensor rotates relative to the orientation of the soil fabric. The observation and experiment on flow-type behavior of liquefied soil has revealed that the influence of inherent fabric anisotropy on the residual strength of a granular soil is so drastic that the inherent anisotropy can no longer be ignored in sand modeling.
The influence of fabric anisotropy has been known since the emergence of the Geo-mechanics [6] . Generally, all particulate materials, particularly sand, are naturally anisotropy and formed by a process of sedimentation or water/air bloviated. However, the process of sedimentation, as well as the pressure exerted by the overlaying material, causes flat soil particles to be oriented with their longest dimensions approximately parallel to the plane of deposition. This condition is achieved upon the minimum energy law of nature and causes a better contact of particle with the others in vertical direction such that the soil strength against lateral expansion will be more than the other directions. As a result, similar to the effect of lateral pressure the soil strength increases against applying loads in perpendicular direction to the bedding plane. However, most natural soils possess a transversely isotropy structure. The sedimentation process may also lead to the formation of alternating layers of different texture. Such a stratified arrangement is typical in natural deposited soils. However, anisotropy manifested itself through the directional dependence of the deformation characteristics of granular materials and this has been widely documented in the literature [7,8] . Furthermore, the degree of anisotropy may vary quite significantly, depending on the soil composition, or even other sources such as electro-chemical properties of pure water, consolidation history, etc.
Given the intrinsic oriented nature of soil fabric, it is important to include the effect of anisotropy in a rational way. The major obstacles, in this respect, or our ability to properly define the spatial and temporal variations of the soil properties, deformablities, hardening and boundary conditions. The value of a model lies primarily in its ability to capture the basic trends in the material behavior and thereby provide a more realistic representation of the problem.
However further than the fabric property of natural soil, the response of granular materials to a given stress depends on the orientation of that stress, nevertheless, the particle alignment constituted in rivers, beach, coastal sands or artificially deposited sands. This nature is also observed in a random arrangement of constituent glass balls which obviously there is not the effects of fabric on the composition by Oda, 1977. This feature is in general known to be due to gravity and the motion tendency of the constituent particles towards the earth center. This tendency can actually affect the geometry of contact points of each particle no mater of shape and particle sizes.
So far, many investigators tried to carry out reasonable relations for the effect of other parameters on internal friction angle [9,10] . The conclusion of their investigations can be expressed as follows. The internal friction angle is of particular importance, whatever its physical meaning; it is not constant for a given soil, even though the void ratio or relative density is fixed, but rather changes greatly depending on three other factors. The first factor is intermediate principal stress σ ' 2 or alternatively b=(σ ' 2 -σ ' 3 )/(σ ' 1 -σ ' 3 ), i.e. b=0 in a triaxial compression test (σ ' 2 = σ ' 3 ), b=1 in a triaxial extension test (σ ' 2 = σ ' 1 ) and b=0.2-0.4 in a plain stress test. The second factor is the effective confining pressure P ' c or effective mean stress (σ ' m ) and the third factor is the principal stress axes direction with respect to the common trend that φ′ increases rapidly, by 10 to 20% in the case of dense specimens, as b increases from 0 (triaxial compression ) to 0.2 to 0.4 (plane strain) and beyond the b value for plain strain, φ′ decreases gradually.
On the other hand, since granular materials (sand, gravel or even glass beads) show more or less anisotropy in their micro-fabric, their internal friction angles commonly change depending on the major principal stress axes direction with respect to the fabric anisotropy [7,11] . The inherent anisotropy is produced by the tendency for non-spherical particles, when deposited under the action of gravity, to lie in the horizontal plane (called the bedding plane) with the long axes parallel to it. This type of anisotropy is well preserved until failure, so it has an appreciable effect on the local spatial distribution of particle rearrangement as well as activated internal friction angle.
In other words, it can be said that this kind of anisotropy diminishes due to particle rotation and movement during shear deformation around a certain deformed location. It must be noted here that this concept is different from the concept of induced anisotropy. Induced anisotropy is generally, initiated and constructed during plastic shear deformation and plays a key role in understanding the plastic behavior of granular soil in a general stress state, including the rotation of principal stress axes [12] .
Another kind of anisotropy is the one that initiated by the pore water pressure during undrained tests. The main source of this anisotropy is known as the variation of the pore water pressure coefficient of. This anisotropy is primarily responsible for the observed difference in the undrained strengths of horizontal and vertical specimens, depending only slightly, if any, on the anisotropy of the shear strength parameters in effective stresses.
In this study, a certain function of local strength variation is defined that depends to overall causes of inherent anisotropy at that location. This function can present the maximum shear to normal stress as tan(ϕ) in any probable sliding direction through the medium. Consequently, any change in major principal stress axis direction is faced with new sets of maximum strengths against the sliding mechanism of grains.

Strain distribution around a point:
In general continuum mechanics, to define the strain distribution at a point, the components simply are considered on the outer surface of a typical dx, dy, dz element. This method makes the solution to consider uniform and homogeneous strain distribution of components over the outer surface of such dx, dy, dz element on three perpendicular coordinate axes. There is a further consideration in addition to the requirement that displacements of a granular medium provided due to slippage/widening/closing between particles that make a contribution to the strain in addition to that from the compression of the particles. Consider two neighboring points on either side of the point of contact of two particles. These two points do not in general remain close to each other but describe complex trajectories. Fictitious average points belonging to the fictitious continuous medium can be defined which remain adjacent so as to define a strain tensor. The problem presents itself differently for disordered particles compared with the orders sphere of equal sizes. In this case small zones even may appear in which there is no relative movement of particles. This can lead to specific behavior such as periodic instabilities known as a slipstick that creates non-homogeneity in strains and displacements.
The effects of non-homogeneity in mechanical behavior of non-linear materials are very important and must somehow be considered. Furthermore, these nonhomogeneities mostly are neglected in mechanical testing because strains and stresses are usually measured at the boundary of the samples and therefore have to be considered reasonably within whole volume.
Solving a nonlinear problem, the mechanical behavior depends strongly on the stress / strain path as well as their histories. Upon these conditions, it may be claimed that the consideration of strain components along three perpendicular coordinate axes may not reflect the real historical changes during the loading procedure. In the most extreme case, the definition of a sphere shape element dr (instead of dx, dy, dz cube) carrying distributed strain similarly on its surface can reflect strain components on infinite orientation at a point when dr tends to zero.
The finite strain at any point in three dimensions by coordinates (x, y, z) relates to the displacements of the sides of an initial rectangular-coordinate box with sides of length dx, dy and dz to form the sides of a parallelepiped as shown in Fig. 1a. This configuration of strain is established by considering the displacements of the corner points (x, 0, 0), (0, y, 0) and (0, 0, z). This kind of strain approach leads to define a (3×3) strain tensor including six components to present displacement gradient matrix at a node. Accordingly, any displacement and corresponding gradient have to be defined as independent components on three perpendicular coordinate axes. Figure 1 shows two rectangular and sphere elements and a typical deformed shape of them. Obviously, there is a certain history of displacement on any random orientation through the element. These are abbreviated in three when a box -shape element is employed. To avoid not missing any directional information of strain, a spherical element carrying strain components over its surface, as tangent and normal to the surface must be employed. This form of strain certainly represents a better distribution includes all directional information. Certainly, to obtain the strain components as presented on planes around box element, strain variation is integrated over the sphere surface. However, a predefined numerical integration may be employed to ease the solution. Numerical integration generally simulates the smooth curved sphere surface to a composition of flat tangential planes make an approximated polygon to sphere surface. Any higher number of sampling plans, the approximated surface is closer to the sphere. Clearly, if the number of sampling plans has taken six, the approximated surface is the same as normal dx, dy, dz box element.
Multi-plane framework: Granular materials consisting of grains in contact and surrounding voids are particulate media that mostly are considered a continuum for ease. The accurate behavior of such particulate materials is to be investigated through micro-mechanics. However, the micro-mechanical behavior of granular materials is therefore inherently discontinuous and heterogeneous. The macroscopic as an overall or averaged behavior of granular materials is determined not only how discrete grains are arranged through the medium, but also by what kinds of interactions are operating among them. To investigate the micro-mechanical behavior of granular materials, certainly, the spatial distribution of contact points and orientation of grains must be identified. At the engineering point of view, the main goal is to formulate macro-behavior of granular materials in terms of micro-quantities. However, there exist two well-known theories that explain the relation between micro-fields and macro-fields as macro-micro relations, in a consistent manner as the average field theory and the homogenization theory.
For a granular mass such as sand that supports the overall applied loads through contact friction, the overall mechanical response ideally may be described on the basis of micro-mechanical behavior of grains interconnections. Naturally, this requires the description of overall stress, characterization of fabric, representation of kinematics, development of local rate constitutive relations and evaluation of the overall differential constitutive relations in terms of the local quantities. Represents the overall stress tensor in terms of micro level stresses and the condition, number and magnitude of contact forces has long been the aim of numerous researchers [13,14] .
Multi-plane framework by defining the small continuum structural units as an assemblage of particles and voids that fill infinite spaces between the sampling plans, has appropriately justified the contribution of interconnection forces in overall macro-mechanics. Upon these assumptions, plastic deformations are to occur due to sliding, separation/closing of the boundaries and elastic deformations are the overall responses of structural unit bodies. Therefore, the overall deformation of any small part of the medium is composed of total elastic response and an appropriate summation of sliding, separation/closing phenomenon under the current effective normal and shear stresses on sampling plans. These assumptions adopt overall sliding, separation/closing of inter-granular points of grains included in one structural unit are summed up and contributed as the result of sliding, separation/closing surrounding boundary planes. This simply implies yielding/failure or even ill-conditioning and bifurcation response to be possible over any of the randomly oriented sampling plans. Consequently, plasticity control such as yielding should be checked at each of the planes and those of the planes that are sliding will contribute to plastic deformation. Therefore, the granular material mass has an infinite number of yield functions usually one for each of the planes in the physical space. Figure 2 shows the arrangement of artificially polyhedron simulated by real soil grains. The created polyhedrons are roughly by 13 sliding planes, passing through each point in the medium. The location of tip heads of normal to the planes defining corresponding direction cosines are shown on the surface of unit radius sphere.
In the ideal case, the normal integration is considered as summing up the individual micro effects correspond to an infinite number of micro sampling plans. The choice of 13 planes for the solution of any three dimensional problems is a fair number. The orientation of the sampling planes and direction cosines of two perpendicular on plane coordinate axes and weighted coefficients [15] for employed numerical integration rule and calculation of the stress tensor of each plan are shown in Fig. 3.
Upon yield criterion in plasticity, the stress condition exceeds the yield limits, plastic sliding or widening/closing take place as an active plane. Therefore, one of the important features of the multi-plane framework is that it enables identification of the active planes as a matter of routine. The application of any stress path is accompanied with the activities of some of the 13 defined planes at any point in the medium. The values of plastic strain on all the active planes are not necessarily the same. Some of these planes initiate plastic deformations earlier than the others. These priorities and certain active planes can change due to any change of direction of stress path, a number of active planes may stop activity and some inactive ones become active and some planes may take over others with respect to the value of plastic shear strain. Thus the framework is able to predict the mechanism of failure.
Strength ellipsoid: In general, quantitative description of initial micro-fabric would enhance characterization and forecasting of sand behavior under different loading. On loading the fabric is continuously altered. Hence, it is necessary to develop techniques to quantify the changes in fabric as well. While the material distorts the fabric of the material changes and so there is the strain or a displacement field in the material. Consequently, the strain and induced fabric of a material are inherently related to each other.
A popular approach for formulation of strength criteria for anisotropic granular materials is the generalization of isotropic ones. Such a criterion is usually geometrically interpreted as a limiting envelope in a stress space; which means that a condition of failure occurs when a given stress vector touches the failure envelope. Since the condition for failure is intrinsic to the material, the failure criterion can be defined differently for any probable sliding plane through the material. Accordingly, the stress ratio cannot exceed the corresponding value of (tan (φ)) on neither on the planes of weakness nor on any other plane which not tending to slide.
On a loading orientation inclined by three angles (direction cosines) respect to the bedding plane, a certain sliding mechanism composed of active sliding planes provides a value of stress ratio which corresponds to the most active plane and has a limitation of (tan(φ)) that is governing sand strength against sliding. Therefore, on any orientation within the sand, the state of strength depends on the geometry of bedding plane and the orientation of applied load with respect to the bedding plane. To describe the strength, (tan(φ )), at any orientation, it is needed to find a way of summarizing the configuration of different strengths corresponding to all of probable directions passing through the medium.
For ideal granular media with no preferred orientation a spherical envelope of strengths (tan(φ)) may provide uniform sliding strength on any orientation. However, to consider fabric effects due to bedding plane, an ellipsoidal envelope of strength may be the most suitable presentation of strength variation on different directions. The longest diameter of this ellipsoid has been always normal to the bedding plane. Configuring the 13 predefined planes in strength ellipsoid provides a certain elliptical section on each plane that present the variety of strength with respect to sliding orientation. In other words, the tips of the arrowhead of strength value of different orientations, collectively define a built up geometrical surface called the strength ellipsoid. Obviously, the size of strength ellipsoid of each plane is different and present maximum and minimum strengths for sliding along the longest and shortest ellipse respectively. The other sliding orientations are faced with strength limitations in between depend on the direction of shear stress on plane respect to the bedding plane.
Adopting the multi-plane mechanism of sliding planes configured in Fig. 3 with respect to orientation of applied major principal stress axis, these planes are configured in symmetric manner around major principal axis. Any change in principal stress axis direction creates a new set of strength ellipses with different strength against sliding directions on different planes.
Extensive similar aspects to this subject were proposed by different authors [16][17][18][19] . Accordingly, the emphasis of these studies, we're on the evolution from inherent to induce anisotropy. Oda [7] proposed a diagram of frequency of contacts as a function of the orientation of their normal in the direct shear test, before and after shearing as fabric ellipsoid.
This study carried out for micro-fabric behavior of loose, medium and dense granular materials led to the establishment of a statistical criterion of natural anisotropy based on hypotheses that experience accepted as probable. The fabric anisotropy law is represented as a spatial closed ellipsoid strength function in x, y and z coordinates as follows: 2 2 2 x / B y / B z / A 1 0 A, B and C are three mutually perpendicular diameters of the ellipsoid respectively. A construction of typical ellipsoid is shown in Fig. 4.
Furthermore, to overcome the anisotropy concerns the loading orientation, a possibly having all different probable sliding mechanisms must be provided in the used model. In this way, the application of any arbitrary loading or stress path leads to a certain sliding mechanism obeys the minimum energy level in natural law. These possibilities are provided with an elasticplastic Multi-Plane model.
To find the strength ellipsoid diameters, two triaxial standard compression tests must be arranged; one with horizontal bedding plane and vertical loading axis (test one) and the other with both vertical bedding plane and loading axis (test two). According to author experience, in both of two test planes 1 to 4 are mostly active, however, in the second case because of 90 degree rotation of load on axis, the smaller ellipse of strength on active planes are provided by the main strength ellipsoid. Their less friction angles are applicable while loading axis diverged respect to the normal axis to the bedding plane. In triaxial compression tests, the sliding orientation of active planes is along the longest diameter of each plane strength ellipse. Furthermore, these four most activated planes in test one are symmetrically located around loading axis that is identical with normal line to the bedding plane. Therefore, for any axe-symmetry loading conditions similar to triaxial compression test, there should be axesymmetry of strength respect to the normal axis to the bedding plane. In this case, the obtained ϕ value of tests and the coordinates of tip head arrows of shear strain on planes provide a unique equation relates the unknown strength ellipsoid parameters. However, due to axisymmetry, two minor diameters of strength ellipsoid are equal, so, B=C. To find the coordinates of tip head arrows of shear strain, one simply can obtain direction cosines of shear stress/strain on corresponding plane, considering the length of the arrow to be equal to tan(ϕ), one point of strength ellipsoid is known. It must be added that in triaxial standard tests, the coaxiality of strain and stress is valid, therefore, to obtain any active plane sliding that is enough to find the direction of corresponding shear stress. In the second compression test, the loading axis is rotated by 90 degrees, so, the same active planes rotated by same angle and their strength ellipses as well. In this case, the geometry of strength ellipsoid is the same as test one. However, different tan(ϕ) and coordinates of tip head arrows are provided that lead to provide the second relation between strength ellipsoid parameters. The simultaneous solution of both equations presents the unknown parameters A, B and C. Assuming the direction cosines of the advanced active planes in first and second tests as , , m n ℓ and ' ' ' , , m n ℓ , respectively, A and B are calculated as follows: The simultaneous solution of two equation yields: n tan ( ) n tan ( ) Performing a compression plane stress, axesymmetry condition is not available any more. Prevention of out of plane strain affects the activation and sliding orientation of planes. In this case, the conservation of the minimum level of energy law forces the mechanism to occur in a different manner as well as geometry. The change in sliding orientation on an active plane from the first natural possible case may make a necessity of grain rolling under constrained conditions. This may lead the sliding to face on local interlocking that generally appears as unstable higher friction angle ϕ f that creates an unstable larger strength ellipsoid for these kinds of sliding planes. This unstable ellipsoid will disappear as soon as the stress path passing over top peak shear strength and softening mode and strength condition comes back to normal strength ellipsoid. Certainly, the larger strength ellipsoid provides a larger unstable intersection strength ellipse on corresponding plane.
To find out the values of internal friction of 13 planes oriented inside a certain strength ellipsoid, first the direction cosines of stress vector as ' ' ' , , i i i m n ℓ are calculated as follows: The value of ' tan( ) The direction of calculated shear stress on i th plane is associated with a certain value of internal friction angle ϕ i in strength ellipsoid. This friction angle can be obtained through the equation of an intersected ellipse plane with strength ellipsoid having the direction of shear stress direction cosines. Simply, any change of shear stresses on planes, is faced on new sliding mechanism and strengths.

Constitutive equation:
Constitutive modeling of particulate material including different features has been the subject of numerous investigations during the recent years, primarily because of the increasing awareness of the complexity of the loading conditions to which particulate material structures are subjected and the corresponding need for more accurate analysis for prediction of safety of such structures. The parallel development of more powerful and efficient numerical methods of analysis has motivated and allowed the use of sophisticated constitutive models beyond the linear or simple non-linear elastic-plastic constitutive laws which were utilized in the early stages.
Most models proposed are based on the theory of elastic-plasticity, incorporating different yield criteria, flow and hardening rules. Strain hardening models according to various isotropic, kinematics or mixed hardening rules have been proposed. These models usually deal with a single or a combination of stress invariant. Principal axes rotation of stress/strain or both and induced/inherent anisotropy have been observed in many tests. However, a model based on invariant of stress/strain tensors, therefore cannot cope with the real behaviour of particulate material such as soil under a complex loading program while either the values of stress or strain invariant are kept constant.
The first multi-laminate model presented by Zienkiewicz et al. 1977. A multilaminate model for granular material was developed by Sadrnejad , et al. 1987 andSadrnejad, 1990. Also, a micro-plane model was developed by Bazant, et al, 1983. For the soil mass, the overall stress-strain increments relation, to obtain elastic-plastic strain increments (dε ep ), is expressed as: C ep is elastic-plastic compliance matrix. Clearly, it is expected that all effects and changes in elastic and plastic behaviour be included in C ep . To find out C ep , while decomposing the global behaviour into different sampling planes, the constitutive equations for a typical slip plane must be considered in calculations. Consequently, the appropriate summation of all provided compliance matrices corresponding to considered slip planes yields overall C ep , therefore, strain increment at each stress increment is calculated as follows: Lε and L σ are transformation matrices for strain and stresses, respectively.

Constitutive equations for sampling points:
A sampling plan is defined as a boundary surface that is a contacting surface between two structural units of polyhedral blocks. These structural units are parts of a heterogeneous continuum, for simplicity defined as a full homogeneous and isotropic material. Therefore, all heterogeneities behaviour is supposed to appear in inelastic behaviour of corresponding slip planes. In many cases, however, the medium is known to be heterogeneous and the notion of a continuum is used to describe it on a scale very much larger than the scale of the real particles. When this approach of 'presmoothing' is taken a priori, without any knowledge of the distribution and aggregation of specific microstructure, information on all internal detail, on the distribution of inter-granular stresses, strains and many other real features is forfeited. Since in reality, this information is necessary to understand the overall deformation resistance of the soil, this aspect becomes too complicated. Therefore, the material that is contained inside a structural unit is treated as a 'black box'. In many instances, the scale of the microstructure is coarse enough to be out of the range of such specific considerations of slip theory and the individual component blocks can be considered as a continuum with well-defined plastic resistances and hardening behaviour. In this research, the individual component blocks of the overall media deform collectively as a heterogeneous (but compatible in deformations with other blocks) assembly of continua, interacting with each other only through the boundary conditions applicable at their various interfaces. The deformations of such coarse heterogeneous assemblies are best considered in full detail, preserving the information on the internal variety of effective deformational resistances in individual component blocks and associated internal stresses. This can then be followed by averaging or 'post-smoothing' approach permits the monitoring of the evaluation of internal deformations in addition to the over all deformation resistances.
The elastic-plastic constitutive equations of a sampling plane start with the classical decomposition of strain increments under the concept of elastic-plasticity in elastic and plastic parts are schematically written as follows: The increment of elastic strain (dε e ) is related to the increments of effective stress (dσ ) by: [D G io and h are material properties, e is the void ratio and Pa is the atmospheric pressure. The bulk modulus may be computed by assuming a value for Poisson's ratio. The dilatancy factor for i th plane is defined as: C′ i is cohesion of soil for i th plane. A loading parameter, L i is defined for i th plane as follows [20] : is calculated as follows: Therefore, the global compliance matrix based on multi-plane numerical integration is found as follows:

Verification of model by toyoura sand:
The test result on Toyoura sand, whose specific gravity is 2.66 g/cm 3 , was employed for verification of proposed model. The maximum void ratio e max of this sand was 0.977 and e min was 0.597. The grain size distribution curve is shown in Fig. 6. The mean grain size D 50 was 0.20 mm and the uniformity coefficient was 1.24. Linear variations of e ph and e f versus (P/P a ) ξ , as presented by Li et al. [21] are adopted.

IDENTIFICATION OF PARAMETERS AND RESULTS
For an axe-symmetry anisotropy of soil that is more practical, two triaxial standard compression tests upon horizontal and vertical bedding plane is enough to identify A, B and C. The strength and dilation parameters of multi-plane model are G 0 , ν, φ fi , η ci , e Γ , C λ , ξ , d 0 , h 1 , h 2 and n. The first two parameters correspond to the elastic behavior of soil skeleton. φ fi ,is obtained from the built up strength ellipsoid on corresponding bedding plane for 13 planes as presented in Sadrnejad [6] . The ratio of η ci /tan(ϕ i ) is assumed initially the same for all directions. This ratio is obtained during the calibration of stress-strain of triaxial standard tests. This parameter assumed to vary after exceeding phase transformation line to exceed one. This variation means that the strength ellipsoid approaches to sphere at failure condition. Certainly, there should be a continuous reform of this ellipsoid due to the rotation of principal effective stress axes, however, to ease the solution, this variation can be neglected up to the steady state line. Accordingly, this reform has been considered during post-liquefaction. Starting from the hydrostatic initial stress condition and no plastic strain history, these parameters except friction angles are assumed to be the same for all planes and are to find by numerical calibration with the results of the stated two triaxial compression tests.  (tan(φ f )) max =0.8665 and (tan(φ f )) min =0.7602. A=0.8563, B=C=0.7752 [6] .
To present the ability of proposed model, the test results conducted by Verdugo & Ishihara on Toyoura sand are produced by model. The obtained model results through the use of calculated parameters are presented as the comparison of stress deviator versus effective mean stress (stress paths) and stress deviator versus axial strain. Figure 7a and b, shows the comparison stress paths of model result with tests for dense Toyoura sand as stress deviator versus different initial mean stress. Also, the comparison of the stress deviator versus axial strains are shown in Fig. 8a and b. The comparison of the stress path to loose Toyoura sand are shown in Fig.  9a and b. The comparison of the stress deviator versus axial strain for loose sand are shown in Fig. 10a and b.
To show the capability of proposed model in predicting the effects of bedding plane and canges in orientation of principal stress axes, the hollow cylindrical test results presented by Nakata, et.al. [22] were produced and compared with test results. Figure  11a to 11f are showing like to like comparisons.
These comparisons show that the model has predicted tally results as the effects of principal stress axes orientation on mechanical behaviour of Toyoura sand.

CONCLUSION
A multi-plane model incorporating the steady/critical state concept was developed for the undrained behaviour of sand. The feature of this model is the use of a pressure dependent peak stress ratio parameter that approaches the steady/critical state value as the steady/critical state normal stress is approached. The model showed to be capable of accurately predicting the undrained behaviour of sand over a wide stress region.
A method to solve anisotropy of soil as the effects of natural micro-fabric and also due to inclination of direction of applied load with respect to bedding plane is presented and verified. This method is simple and based on the minimum energy level in natural law. Despite of the direction effect on soil strength, the presented distribution of strength by rotation of bedding plane is of a unique form through the material. This aspect simplifies the use of the presented method to find strength ellipsoid parameters.
A micro-plane numerical algorithm is also presented for a better anticipation of load inclination effects through material. In this way, the directional information and effects of applied load orientation on mechanical behaviour of material are addresses and considered. Further than possibility of predicting inherent anisotropy, this rational way facilitates to the model to predict the effects of stress/strain principal axes rotation, induced anisotropy and a potential to solve anisotropy of material through defining different mechanical behaviour on different orientation. This is achieved by the use of a generally simplified, applicable, effective and easily understandable relations between micro and macro scales. These relations demonstrate an easy way to handle any heterogeneous material property as well as the mechanical behavior of materials.
This model is able to solve a three dimensions plasticity problem by a rather simple theory based on the phenomenological description of two dimensions plastic deformation and kinematics hardening of materials. This, is actually, achieved in such a way that the application of some difficult tasks such as induced and inherent anisotropy and rotation of principal stress and strain axes where there may be not co-axiality of them during plastic flow, to be predictable. Accordingly, the sampling plane constitutive formulations provide convenient means to classify loading events generate history rules and formulate independent evolution rules for local variables.