Simulation of shear bands with Soft PARticle Code (SPARC) and FE

The aim of this paper is to numerically investigate the development, thickness and orientation of shear bands, in biaxial test with two approaches towards solving problems of continuum mechanics, namely the meshless “Soft PARticle” method and the mesh based Finite Element method. Soft PArticle Code (SPARC) is a straightforward collocation numerical method based on strong formulation, in which a first order polynomial basis is adopted for the evaluation of spatial derivatives in partial differential equations. A novel nonlinear constitutive model— barodesy for clay, is adopted in this study. The biaxial test, which involves homogeneous, and later inhomogeneous localized deformation is simulated using the Soft PArticle Code and the Finite Element method. The inclination and thickness of the shear bands are evaluated and analysed with the earlier experimental, theoretical and numerical investigations. Furthermore, simulation results are compared and presented to demonstrate the advantages and limitations of SPARC in comparison to FE method.


Introduction
Meshless methods have drawn attention in the past decades due to their simplicity in numerical formulation and ability of modeling large deformations. Compared to mesh B S. M. Iman Bathaeian iman.bathaeian@uibk.ac.at; iman.bathaeian@gmail.com 1 Division of Geotechnical and Tunnel Engineering, University of Innsbruck, Technikerstr. 13, 6020 Innsbruck, Austria based methods such as FE methods, a true meshless method maintains neither the connectivity between the nodes nor uses any background mesh, therefore, no mesh is required. In addition, the calculation of spatial derivatives in partial differential equations is achieved by means of the information (the parameter of interest e.g. stress or velocity and the position vector) carried on by each particle. The point collocation method (Aluru 2000) and finite point methods (Onate et al. 2001), for instance, belong to this category. These methods use a finite support domain (small in comparison to the whole study domain) to approximate the spatial derivatives, and hence belong to collocation methods, which result in significant reduction of computational efforts (Zhang et al. 2000). Thereafter, the governing equations are solved, meaning that such a formulation is "strong".
It has been shown that although SPARC adopts basic polynomial linear approximation method, it is capable of modeling strain localization (Chen 2005). In addition, the formulation for stress boundary condition is straightforward.
Finite element (FE) methods have been developed since the 1950's. They are often regarded as benchmarks for new methods (Liu et al. 2006b, a) due to their robustness and accuracy. FE uses the weak formulation to solve the underlying differential equations. In this paper, the simulation results of a biaxial test using SPARC are compared with those obtained with an FE method. Two sets of biaxial test on overconsolidated specimen, once with a weak zone and once homogeneous are modeled. The weak zone guarantees the occurrence of inhomogeneous deformation and localization of deformation after the homogeneous deformation and manifests strain localization. Strain localization in real samples has e.g. been captured by means of x-ray on test samples (Kolymbas 1972). Modeling of strain localization is a challenging task, it can thus be used to test the capability of numerical simulations and serves as a good numerical example for comparison herein.
In our simulations, barodesy for clay (Medicus and Fellin 2015) is used as constitutive model. Barodetic constitutive models (Kolymbas 2012a, b;Medicus and Fellin 2015) take into account the influence of stress state T and void ratio e of the material. Their incremental stiffness (Kolymbas 1999) imposes convergence difficulties (Fellin and Ostermann 2002) especially at the occurrence of strain localization. This paper is organized as follows: in Sect. 2, the framework of SPARC is briefly explained, Sect. 3 provides some information regarding the adopted constitutive model barodesy for clay, in Sect. 4 the conventional soil parameters as defined by Mohr-Coulomb are obtained, Sect. 5 presents the simulation setups and results using SPARC and FE method and in Sect. 6 the acquired shear bands are investigated and compared. Conclusions are given in Sect. 7.

SPARC framework
Soft PArticle Code is a straightforward framework for numerical simulation based on collocation strong formulation (Chen 2005;Ostermann et al. 2013). Continuum is presented using particles that carry physical information such as position vector x, velocity vector v, Cauchy stress tensor T, void ratio e and etc. The equilibrium equation for each particle reads: div t+Δt T + ρg = 0 (1) Since in this study the body force induced by gravitational acceleration g is negligible, g is set to zero in our simulations. The velocities of particles are the unknowns in SPARC. An initial guess of the velocity field v k=0 is considered as the start point of the solution finding process (see Fig. 1). The velocity gradient L = grad t v k=0 is then evaluated by means of "first order polynomial" interpolation. We then obtain the stretching tensor D and the spin tensor W, which are the symmetric and anti-symmetric part of L as shown in Fig. 1. The objective rate of the stress tensorT t at time-step t is evaluated by the constitutive model "barodesy for clay",T which is a function of the current stress state T t , e t and D of the particle under evaluation. Due to highly nonlinear nature of barodesy for clay, the fourth-order Runge-Kutta integration is applied to obtain a more accurate objective stress rate T t [see Chen (2005)]. Since the void ratio e is also one of the variables affecting the stress rate, the void ratio is also updated in every sub-step of the Runge-Kutta fourth-order by calculating the rate of void ratiȯ where tr D is the sum of the diagonal components of the stretching tensor. Application of Jaumann-Zaremba objective stress rate yields the stress rateṪ, and finally the stress state and void ratio are updated as follows: The above mentioned algorithm is applied to all particles and the updated stress tensor T t+Δt and void ratio e t+Δt for all particles are obtained. The equilibrium equation Eq.
(1) can be built once the T t+Δt has been calculated. The non-linear system of equations consisting of Eq. (1) is solved by Newton method in order to determine the unknown velocities v t+Δt k+1 , where k is the iteration of the Newton solver until the solution which satisfies the Eq. (1) is obtained. Once the solution v t+Δt is available, the position of the particles are updated, and the solution of the velocity field from the previous time step is used as the first guess for the next velocity field time step.
Since SPARC is a collocation method, the Jacobian is sparse. In case of the presence of stress boundary conditions, such as the hydrostatic pressure applied on the biaxial test specimen, the kinematic boundary condition T t+Δt n t+Δt − n t+Δt p = 0 is applied, where n and p are the surface normal vector and pressure, respectively. Replacing Eq. (1) as the governing equation for boundary particles subject to surface pressure p. It can be seen that the spatial derivatives of the velocity field and stress field (i.e. gradv and divT) and the unit vector normal to the free surface must be evaluated to obtain the partial differential equations (PDE) in Eq. (1) or in Eq. (8). To achieve this, the linear polynomial approximation method is adopted. The following first order polynomial for m-dimensional problems is used: wheref is an approximation to the parameter of interest f . Its coefficients a j are computed in a least-squares sense, benefiting from the information f carried by the adjacent particles, the so-called neighbors. Here comes a 2-dimensional example for the evaluation of the coefficients a j for a particle with index i: . . .
where n n denotes the total number of neighbors of particles i; the subscript denotes the corresponding physical property belonging to particle i; subscripts 1, 2, dots and n n denote the index of neighbors of particles i; F (i) and X (i) are matrices consisting of the parameter of interest f and the position vectors x stored at neighbors of particle i, respectively. The coefficients can thus be obtained with Herein, neighbors of the under-evaluation particle, are particles, which lie within a cut-off radius r of that particle. Note that the particle under evaluation is also counted as its own neighbor. The order of the neighbors in Eq. (11) is irrelevant, but the order of Eq. (10) must correspond to Eq. (11). Once the coefficients off are obtained, the partial derivatives off are also obtained: ∂f In SPARC, the logarithmic definition of strain is implemented so as to account for large deformations. It is believed that although the homogenized deformation of the whole specimen can be considered small enough to fit to the description of infinitesimal strain theory, the shear deformation in the shear bands is large enough to invalidate the theory of infinitesimal strain. The complete algorithm of how SPARC works is summarized in Fig. 1.

Constitutive model
Barodesy is a constitutive model which can be seen as a further development of Hypoplasticity, which does not include notions of elastoplasticity such as yield and plastic potential functions. Barodesy, pioneered by Kolymbas (2012c) was first developed for sand and then developed for clay (Medicus and Fellin 2015). It is formulated in stress-rates (T) as a function of the stretching tensor (D), the actual stress state T and the void ratio e: It is based upon the two rules by Goldscheider (1976), the so-called proportional paths and the asymptotic soil behavior. Deformations with proportional strain paths result in stress paths which approach a proportional stress path. Proportional strain paths are, for example, oedometric or triaxial compression. Furthermore, common concepts of soil mechanics, as e.g. barotropy, pyknotropy or critical state concept are considered Medicus et al. 2016). The complete equation of "barodesy for clay" reads, where R is the function, which links proportional strain paths to proportional stress paths and includes the stress-dilatancy relation. It must be noted that the constants c i are dimensionless. The (R 0 and T 0 ) are the normalized tensors 1 .
where tr D 0 is the sum of diagonal components of the normalized stretching tensor and is a measure of dilatancy.
Functions f and g are scalar functions that incorporate asymptotic states, critical states, the influence of stress level (barotropy) and density (pyknotropy), where σ * is the reference pressure (equal to 1 kPa) and The constants c i are determined by means of soil parameters as follows, ϕ c : critical friction angle N : ordinate intercept of the isotropic normal compression line (NCL) in the ln p − ln(1 + e) plot λ * : is the slope of the NCL κ * : is the slope of unloading line under isotropic compression in ln p − ln(1 + e) plot The calibration of constants for Dresden clay are found in Table 1.
All parameters can be obtained from a consolidated undrained triaxial compression test. Readers are referred to Medicus and Fellin (2015) for more profound details on calibration and structure of "barodesy for clay".

Critical state in barodesy
The R function of barodesy has the major contribution in defining the critical state behavior, all stress directions R 0 for isochoric (tr D = 0) form a fan in the principal stress state (see Fig. 2). This locus coincides practically with the failure criterion of Matsuoka-Nakai for granular materials, which is defined as: where I 1 = tr T, I 2 = (I 2 1 − trT 2 )/2 and I 3 = det T are the first, second and third invariants of the stress tensor, respectively. The material parameter K M N is a function of the critical state friction angle ϕ c : Fellin and Ostermann have proven in Fellin and Ostermann (2013) that the deviation of the locus defined by barodesy from the one defined by the Matsuoka-Nakai criterion is less than 0.12% for a critical state friction angle of less than 40 • . The barodetic model has been implemented in Abaqus as a user subroutine Umat based on the hypoplastic formulation presented in Fellin and Ostermann (2002). The void ratio is an additional state and output variable in the Umat.

Mohr-Coulomb vs. barodesy
The analytical solutions addressing the orientation of shear bands for granular materials, are based on two conventional soil parameters, friction angle ϕ and dilatancy angle ψ (see Sect. 6 for detailed discussion). Element simulation of triaxial test with barodesy for clay for three consolidation pressures ( p ini = 100, 200 and 300 kPa) were conducted and the results are presented in Fig. 3. Where σ 1 and σ 2 are the principal stress components. The acquired values achieved from the element tests were applied

Simulation setup
A biaxial test can be regarded as a plane strain adaption of a triaxial test, to show shear localization in a 2D numerical setup. The deformation is driven by two lubricated loading caps on the top and the bottom of the specimen to allow lateral deformation. The upper cap moves downward compressing the specimen. The specimen is loaded with constant hydrostatic pressure p being applied to the specimen. The resulting boundary conditions are illustrated in Fig. 6a. In our simulations, the initial void ratio e = 0.45 under a cell pressure p = 100 kPa, corresponding to an over-consolidated clay, is adopted. Therefore, a peak in the stress-strain relationship with post peak strain softening, strain localization or formation of shear band(s) in the biaxial test simulation are to be expected. The particle configuration for the study domain in SPARC is shown in Fig. 6b. The surface particles are subjected to cell pressure p. Velocity components v 2 of particles on loading caps are prescribed, whereas v 1 is unknown. Note that in order to prevent the specimen from horizontal translation, the velocity component v 1 = 0 for the particle in the middle-top of the sample is prescribed. The radius r = α · h is used to determine neighbors for all particles in simulations where α has a value of 1.5-1.7 depending on the simulations 2 . The unit vector n normal to the pressure surface must be determined for the computation of Eq. (8). n is also computed using first-order polynomials in a compacted support presented in the previous section. Note that therefore, neighbors of a surface particle consist only of surface particles.
For the FE simulation, the particles shown in Fig. 6 are used for the nodes of the mesh. The boundary conditions are the same in both simulations. For the discretization, 4-node bilinear plane strain elements with four Gaussian integration points were applied.
Two simulation examples, with and without imperfection in the specimen, are presented in the following. The imperfection is implemented by increasing the void ratio of the particles representing a weak zone (Fig. 6b) by 0.02, resulting in relatively looser state and thus lower stiffness. In FE simulations, the weak zone consists of only one element formed by the same nodes. Given a weak zone, the formation of shear bands is expected to initiate from the weak zone (Fellin et al. 2009). In order to investigate the dependency of the shear bands, simulated by SPARC, on the number of particles, three sets of simulations with initially homogeneous setup for number of particles 180, 231 and 299 were conducted. Simulation with an implemented imperfection were done only for 231 particles. All simulations are summarized in Table 2.

Simulations with initially homogeneous fields
Since the initial stress field and void ratio field are homogeneous, the deformation in the sample shall be homogeneous, meaning that the stress strain curve of all particles must be identical and must overlap with the curve obtained from element test result. The constitutive model barodesy for clay is directly used to obtain the stress strain curve by prescribing the deformation matrix  Figs. 7 and 8, respectively. The SPARC simulation results show that all curves overlap with one another and with the element test curve until ε 22 ≈ 4.6% is reached. This implies that the deformation of the sample for ε 22 < 4.6% is homogeneous. Thereafter, the deformation starts to localize at particles, the localization causes numerical error and the continuation of simulation leads to the accumulation of the error. At the beginning of the simulation, |D| has almost the same value all over the specimen. However, the discrepancies in the |D| field can be recognized, even if the deformation field is relatively homogeneous. When the axial strain (ε a ) approaches 4.6%, strains start to localize. In the end, shear bands, revealing a 'v' shape, occur on the bottom of the specimen (Fig. 12c).
The void ratio field after localisation is shown in Fig. 11c, results show that contraction occurs in the whole sample at the beginning. When strain starts to significantly localize, the void ratio in the shear bands exhibit volumetric increases (see Fig. 11c). This trend is expected to occur in a dense granular sample with strain softening behavior.
In Fig. 12a, b, the deviatoric shear strain, is plotted as a measure of localization of deformation. The FE simulation, on the contrary, shows a homogeneous behavior over the whole deformation, see Figs. 11a and 12a. The evolution of the void ration e shows also a homogeneous compression while loading and no localization can be modeled, even at axial strain ε a = 20%.

Simulations with imperfection implemented
The stress strain curves in terms of σ 22 and ε a obtained by SPARC are plotted in Fig. 9. For ε a < 2.0% all curves except for those of particles in the weak zone are in good agreement with those of the element test curve. At ε a ≈ 2.3%, strains start to localize significantly in a shear band initialized by the weak zone. Thereafter, strains occur mainly in the shear band. The initial shear band is followed by some other shear bands (see Figs. 11d, 12d) before the program aborts. At this point the solver cannot find any solution even with an extremely small time-steps Δt < 10 −10 . The void ratio field is shown in Fig. 11. Again, contraction occurs in the whole sample at the beginning. However, once the strains start to localize significantly, the changes in the void ratio can be used to demonstrate dilatancy.
The FE simulation shows comparable results. Deformations are homogeneous until the peak. Strain localization occurs at about ε a ≈ 3%, which is larger than ε a ≈ 2.3% by SPARC and closer to the peak of the element test result (Fig. 10). The FE could reach a value of ε a = 20% which is the standard value of the axial strain in triaxial tests. Shear bands start to occur in FE right before the peak of the stress-strain curve and almost in "X" pattern, however, this pattern does not last long and after reaching an ε a of greater than almost 4% the shear strain localizes in a 6 Orientation and thickness of the shear bands Vermeer (1990) conducted theoretical and experimental investigations on the orientation θ and thickness of shear bands in biaxial tests, his investigations show that for fine sands, the orientation of shear bands coincides almost the Mohr-Coulomb solution θ C = 45 • + ϕ/2 and for coarse sands, the Roscoe solution of θ R = 45 + ψ/2 is observed. Where ϕ and ψ are the friction and the dilatancy angles, respectively. Investigations of Han and Drescher (1993)  magnitude of the confining pressure. As mentioned in Han and Drescher (1993), the shear band inclination angle with respect to the minor principal stress decreases when the confining pressure increases, however, the shear strains increase. Experimental results of Han and Drescher (1993), have shown that at higher confining pressures (almost 400 kPa), the orientation of shear bands correspond better to the solution of Roscoe and the shear band inclination is in general much lower than the one predicted by Mohr-Coulomb. In our simulations, the acquired inclination angle with SPARC are about 39.8 • for the test with initially homogeneous sample and about 45.8 • for the test with implemented imperfection. As mentioned before, no shear band has been observed in simulation with FE for initially homogeneous sample, while for the sample with initial imperfection, the acquired inclination is about 46.8 • .
As discussed in Sect. 4, a friction angle of ϕ = 32.9 • and dilatancy angle of ψ = 7.3 • can be attributed to Dresden clay. Considering the solution of Mohr-Coulomb with θ C = 45 • + ϕ/2, we should be expecting an inclination of θ C = 61.4 • which is a clear overestimation for both the results of FE and SPARC.
However, the solution of Roscoe with θ R = 45 + ψ/2 would predict a θ R = 48.6 • which seems to be a more realistic estimation both for FE and SPARC.
As for the thickness and inclination of the shear bands acquired by FE methods using a Hypoplastic constitutive model, Tejchman and Wu have shown in Tejchman and Wu (1996)  on the spatial discretization. In a further investigation, Tejchman and Bauer (1996) benefit from the results of an extension of the Hypoplastic model for polar continuum with a characteristic length, the so-called mean grain diameter d 50 . Their results show, the thickness of the shear band is the same for a fine and coarse mesh. In order to realistically simulate the thickness of the shear zone within a polar continuum, the size of the finite element should be smaller than 5 · d 50 . As it can be seen in Figs. 11b and 12b, the shear band acquired for FE calculations show an inclination of almost 52 • and a width of 0.028 m which is almost equal to the width of 4 elements.
In Figs. 11c, e, f and 12c, e, f, results of SPARC for a homogeneous setup and for different number of particles are presented. As can bee seen, for the lower number of 180 particles, shear bands are not clearly formed and deformation seems to localize on the two corners of up-left and down-right. This phenomenon can be so explained, that in case of homogeneous setup, the shear bands appear as a result of the accumulated error in each time-step, and with less number of particles, the accumulation stays smaller which can lead to later occurrence of shear bands or no meaningful occurrence of shear bands.
For more number of 231 particles, the shear band has a "v" shape at the middle bottom of the specimen and has a thickness, containing almost 6 particles (see Fig. 12c), while for more number of 299 particles (see Figs. 11e, 12e), the shear band is not as thick as by 231 particles and contains almost 4 particles. Furthermore, for 299 particles, shear bands have a symmetric shape not only along x 1 axis, but also along x 2 axis.
As for the inclination, the acquired shear bands for more particles (299) have a slightly larger angle as those acquired for 231 particles which fits better to the solution of Roscoe.
From the acquired symmetric shape, and the larger inclination angle, it can be deduced that for SPARC, the same as for FE methods, the denser the particles, the better the shear bands can be simulated.

Concluding remarks
In this paper, we have shown that linear approximation method used in SPARC can model the formation of shear bands, with or without implemented imperfection. The comparison with simulation result using the Finite Element method indicates that both models can simulate shear bands and strain localization when an initially weak zone has been implemented, however, FE is not capable of building any shear band starting from homogeneous conditions. As Abaqus uses the weak formulation, the equilibrium is not fulfilled for every single integration point, but for the whole problem. Starting from a homogeneous field leads to a homogenous solution for the whole field and for all all calculation steps. As a result, no shear bands can develop. SPARC uses the strong formulation to solve the differential equations. Therefore, the equilibrium is fulfilled at every single particle with a prescribed tolerance. SPARC is capable of simulating shear bands even when the specimen has an initially homogeneous setup, this is due to the numerical inaccuracy and error accumulation in the domain. However, this corresponds the reality for experiments with homogeneous setup, that shear bands still occur.
Simulations with different number of particles with SPARC have demonstrated that the density of particles also plays a role in shape, thickness and orientation of shear bands and the denser the particles, the better the shear bands can be reproduced in the framework of SPARC.
Results of simulations with SPARC and FE with implemented imperfection and the same number of particles or nodes, have shown that the inclination and shape of shear band from both methods is comparable and almost similar. One can deduce that simulation of shear bands is not strongly dependent on the applied numerical method, but mainly on the density of particles, mesh size, or the implemented constitutive method, as discussed e.g. in Sect. 6 for constitutive models, developed for polar continuum. The idea of using the velocity gradient and the stretching tensors in the framework of SPARC, provides a more comprehensive structure for dealing with large deformations such as strain localization in granular materials. Furthermore, it offers the advantage that the stretching tensor can be directly used as an input of many constitutive models such as hypoplasticity and barodesy. As for further numerical investigation of shear bands, it is aimed to develop SPARC for saturated conditions and conduct the same biaxial tests for hydromechanically coupled conditions and investigate the development of shear bands for both dense and loose samples and to compare the results with the experimental investigations on undrained biaxial tests.