A Constitutive Model for Multi-Line Simulation of Granular Material Behavior Using Multi-Plane Pattern

Problem statement: Improper understanding of material behavior prevents the efficient and correct usage of available materials and consequently, increases the construction and maintenance costs and even unsuitable construction. Considering the necessity of exact investigation about material behavior, several researches have been carried out in this field but the majority of these researches did not propose a general method for prediction of granular material behavior. Furthermore, many of the methods proposed by researchers are not able to present the properties such as the orientation of failure mechanism of propagating plasticity in materials. Approach: In this study, a general method was proposed for multi-laminate simulation to predict the behavior of materials. The general applicability of this method for prediction of granular material is one of its significant advantages. The study was carried out in the framework of multi-plane pattern which is able to predict anisotropic behavior, consider the effect of stress and strain axis rotations in plasticity, consider the semi-micro mechanical history and finally predict the orientation of failure mechanism. The method was presented in a matter that there is no limitation for different shapes of stress-strain curves. Results: It was concluded that using this method, fundamental properties of material such as material fracture, orientation of failure, anisotropic behavior of material, separation of behavior in several planes and rotation of principle axis of stress and strain during nonlinear behavior can be determined. Conclusion/Recommendations: This method can be used for complicated material behavior simulation under seismic loading, cyclic loading or fatigue effects. For future works, the method can be extended by increasing the number of planes. Higher-order equations can also be used to have a more accurate approximation of stress-strain curve.


INTRODUCTION
Misunderstanding of material behavior and occasional predictions with limited accuracy has caused several problems in analysis and construction in different engineering branches. The use of superfluous safety factors to compensate the lack of knowledge will cause great costs in industry. Therefore, using the methods which are able to predict material behavior will enable the users to recognize all of the affecting factors and properly consider the effect of each factor. Semi-linear or multi-linear analysis is one of the simple methods presented for modeling the material behavior [1] . Koushi and Green elastic models, Hypo elastic models, variable-parameter models are some instances of such model. The notable point in these modeling methods is that in some cases, for better understanding of material behavior, it is necessary to change the parameters in a natural and specific manner and the determination of the parameters may be difficult and even impossible. Besides, the limited ability of these models is in the case of simple stress paths. This fact forced researchers to consider more complicated mathematics with various different components in the new models. But the main problem, which is the dependency of material behavior to stress path and stress history [2] still exist. Because of the mentioned reasons, in this study it was tried to propose a process based on a multi-line method in the form of multi-plane pattern which do not have the mentioned problems. In addition, the significant advances in computer sciences have enabled researchers to entrust the computation of different components to computers in an organized system shape. In this manner, more complicated behavioral parameters can be introduced to models. A significant characteristic in multi-plane pattern is the ability of behavior analysis in specific directions. Existence of seems in stones or position of aggregate positioning in sand and alluvium sedimentation and similar cases is caused the material to have different behavior in different directions [3][4][5][6] . Multi-plane pattern makes it possible to determine the material behavior in different directions.
Multi-plane pattern: The basis of multi-plane method is to determine the numerical relationships between inter-particle behavior (micro behavior) and engineering mechanical characteristics (macro behavior) in the form of an essential equation which is obtained by solution of numerical integration [1,6] . In other words, in this case, behavioral properties of material behavioral properties and stress-strain behavior of soil can be determined by investigation of interparticle behavior. Granular materials are composed of indefinite number of solid particles which are in contact to each other and the reaction between particles is due to the force which exists in contact surface. The analysis of particle behavior and their contact surface depends on number, size, shape, roughness and the strength of particles in these surfaces. Consequently, investigation about these materials under mentioned condition is more complicated than the case which has continuous circumstance. Besides, in a simple view, behavior of granular material can be assumed as the combination of particle elastic behavior and plastic slippage in contact surfaces. In this manner, in artificial case, the three dimensional behavior of granular materials can be determined by considering indefinite sample planes in which slippage is possible. These planes divide the material block into a collection of multilateral parts which are next to each other. When a little shear is applied to a multilateral part, elastic shear deformation occurs. If the shear be increased and achieve a specific limit, multilateral parts move along a boundary surface called slippage plane. With increasing the deformations, necessary shear stress for occurrence of more deformation is increased. The total shear stress at any time is equal to the summation of elastic shear deformation in multilateral parts and shear deformation due to the slippage of adjoining parts. When the stress is decreased, the elastic deformation disappears. When the stress is more decreased and achieved to a specific limit, multilateral parts begin to slip in opposite direction. The required shear stress for slippage depends on normal stress and the slippage only occurs when the stress condition passes the yield stress. Moreover, slippage just occurs in slippage planes with directions shown in Fig. 1.
Numerical concept of multi-plane theory: In this theory, the initial basis is calculation of numerical integration from a specific mathematical function on a sphere of unity radius. This mathematical function can present the changes of physical properties on the surface of sphere. The surface of a sphere can be estimated by indefinite number of flat planes which are tangent on different points of the sphere. In this manner, any of mentioned planes have one contact point with the sphere.
By limiting the number of these planes, the number of contact points or basis points can be adjusted and the calculation of numerical integration gives the amount of the function in the mentioned points. Numerical integration of continuous function f (x, y,z) on the surface of sphere can be calculated as the summation of function value F at different points which are multiplied to corresponding weight coefficients. To decrease the error values, the number of these points should be increased. It can be confirmed that the application of 26 sample points will decrease the error to the sixth order. The below equation shows the relationship between numerical and normal integration: Nodal weight coefficient f (x , y ,z ) The value of F function at the point (x , y ,z ) The position of the 26 points on sphere is shown in Fig. 2. For any point, a plane is defined so that the direction cosines are the same as normal vectors of the plane at contact points.
In this manner, any change on the i plane is concentrically related to it. Since for any of desired points a specific plane is specified, and considering the point symmetry, the 26 points are reduced to 13 points and the surface of half sphere is approximated by 13 planes: The orientations of the 13 planes are shown in the center of cube in Fig. 3. On this basis, if the quality of slippage and opening and closing of each plane is organized, the summation of these slippages, openings and closings compose the method of material movement for a point and by integral summation, total movement or deformation effects can be achieved for a specific point. In this manner, it is required to assume a rule for the value of function i i i, i f (x , y z ) which determines the movements of each plane. Behavior determination of the 13 planes: To determine the behavior condition of each plane and relations between plane behavior and total behavior, the behavior tensor is transferred from global coordinates to each of the 13 planes. Strain vector in global coordinates can be transferred to a three-component vector on each plane consists of normal strain and two shear strains placed on the plane by using transformation matrix.
where, j l is the transformation matrix of jth plane: ll mm nn lm l m mn m n nl n l ll mm nn lm l m mn m n nl n l The parameters n,m,l , n,m,l and n,m,l are respectively the direction cosines of normal to plane direction and two perpendicular direction on the plane. In the Eq. 3, ε is strain vector of six components in general coordinates which is defined as following equation: In Eq. 4, σ is six-component stress vector in general coordinates and c is behavior tensor in general coordinates. Substituting the Eq. 4 in Eq. 3 the following equation is obtained: By this equation, strain vector of planes are obtained. By knowing strain vector of each plane and using basic relation of multi-plane, strain vector in general coordinates which is calculated by summation of planes effect is determined by following equation: where, ' j l is transformation matrix from plane to general coordinates: n 2nn 2nn l lm l m lm l m lm nm n m nm n m nm l.n l n l.n l n l.n Substituting Eq. 5 in Eq. 6, following equation is obtained: In other words, the proportion of each plane in general behavior of plane is as following: The proportion of jth plane in general behavior: In the proposed method, by transferring the behavior tensor into the 13 planes, behavior of granular material is initially determined in each specific direction. In other words, deformations in different directions are ruled. Then by having the rule of deformation in each direction and by numerical integration of multi-plane pattern (based on strain transfer from planes to general coordinates), collective effect of plane behavior on general behavior is obtained. Knowing experimental specimen curves from underground behavior of material and using this method, the model components can be obtained. In other words, due to traffic of strains from general coordinates to planes and vice versa, the components of proposed method are obtained.
A notable point in transferring three-axial experimental curves to the 13 planes is the method of creating stress and deformation on them. As it can be seen in Fig. 3, carried out studies have shown that the effect of the planes 1-4 is to create stress-deformation in the directions normal to plane and both shear directions on the plane. In the planes 5-7, the second shear direction will not become activated. In the plane 8, the shear direction one is inactive. In the planes 9-13, just the plane normal direction is activated and no shear will produced in these planes.
In the above equations, the proportion of each 13 planes in general behavior is determined. If it be necessary, relationships between stress and strain vectors in each plane can be determined: On this basis and using numerical integration of multi-plane pattern, general strain vector is calculated: In this manner, as it was showed in the equation, the proportion of ith plane behavior effect on general behavior is as followings: The proportion of ith plane behavior effect on general behavior: Thus the parametric relation of general behavior tensor and behavior tensor of each plane will be determined.

MATERIALS AND MATHODS
In the following, for a general behavior of material, simulation is carried out using this method to illustrate the usage of the method. Stress-strain carves of this general case in shown in Fig. 4.
Firstly, stiffness and compliance matrices are calculated in each step by equations obtained by presented curves. To prevent the ambiguity of compliance matrix, a very small value is assumed for the steep of lines where the steeps were equal to zero. The aim of writing stiffness matrix was to enable the simple summation of equations based on stress and determination of final compliance matrix in each step. It should be noted that during the solution of problem, the effect of shear stress in producing some proportion of normal stress is considered. In addition normal stress equations are considered both in compression and tension ranges. In the following, the corresponding equations of each step of carve are presented: Part 1:

Part 4:
Having determined stress-strain equations for the 26 planes, the compliance matrix is transferred to general coordinates using transfer vectors. It is obvious that for calculation of transfer matrices, direction cosines should be calculated in normal and shear directions. The values of these cosines are tabulated in Table 1.
It should be noted that for the mentioned problem, since the resultant shear stress is given, direction cosines in one shear direction is just used. Using direction cosines and mentioned equations, transfer matrices of the 13 planes are determined and indicated in Table 2.
Using below equations, general compliance matrix is calculated:  T  T T   T  T T  C  0  0  T  T T  0 C  0  T  T T  0  0 C  T  T T   T  T T   T T  T  T  T  T  T T  T  T  T  T  T T  T  T  T T C  T C T C  T C  T C T C  T T  T  T  T  T  T C  T C T C  T T  T  T  T  T  T C  T C T C  T T  T  T  T  T  T     .
As it can be seen, above equation are very long and the calculation is very time consuming. Thus for simplicity, a program is developed in Excel software which calculates the components of general compliance matrix that are constant for all of the 4 material conditions. Using this program, general compliance matrix for all of the 4 parts of the given problem as discussed above.
To compare the results of the proposed model with real results, the results of an experimental test on a stone specimen in 6 different initial pressures is selected. These results were published in London Geotechnical journal. Considering the 3 carves given for this problem, the curve which is corresponding to initial pressure 100 kPa is the most suitable carve. In the Fig. 5, the results of the test are compared to the result of proposed model.

RESULTS
If the method be properly used and the assumptions be correct, the method can accurately simulate the behavior of complicated materials. The following cases are the advantages of the method: simulation of any material behavior with any stress-strain curve, analysis of fracture, determination of failure orientation, Possibility of application of any failure criteria such as Mohr-Coulomb or Drucker-Prager, possibility of increasing the number of under study planes in multiplane pattern and considering the effects of material behavior in the planes perpendicular to principle planes,

DISCUSSION
The reason of observable differences between the results of the model and the results of experimental test is shown in Fig. 6.
As it can be seen from Fig. 6, to predict the behavior, the line 2 is used instead of curve 1. So in the ith iteration, the values of strain are different by the value i ∆ε . Also in the calculation of bulk strains, since the differences exist in all the three axial strains, the value of error will be accumulated. To prevent unacceptable error value, the number of model lines should be increased. This is shown in Fig. 7. In Fig. 7, the curve 1 is simulated by the three lines as it can be seen and the value of error is significantly decreased.

CONCLUSION
In this study, it is tried to propose a method by which the behavior of granular material can be predicted in a simple and direct form. Using multiplane pattern make it possible to study the behavior of material in different directions. Moreover, the kind of failure can be determined by investigating the stresses in the 13 planes. Furthermore, there will be no limitation for different stress-strain curves using the presented multi-line method and this model has no restriction, since general case is studied and so any complicated and different behavior can be simulated. This method can also be used for study the effect of earthquake and cyclic loading on specimens.