Applications of micropolar SPH in geomechanics

A smoothed particle hydrodynamics code based on micropolar continua for geomaterials is developed for problems involving large deformation and shear strain localization. Two typical geotechnical problems, i.e., biaxial compression test and sand column collapse, are simulated using classical and micropolar model to demonstrate the performance of the newly proposed method. A parameter study is given on the scale effect in the micropolar continua.


Introduction
Smoothed particle hydrodynamics (SPH) was originally developed for simulations in astrophysics [22,37]. After decades' development, SPH has been applied to a vast range of problems in CFD and solid mechanics such as free surface flow [44], muti-phase flows [42], dynamic response [30] and explosion simulations [34]. As a purely Lagrangian meshfree method, SPH is suitable for problems involving large deformation and post-failure, which are challenging for traditional numerical methods like Finite Element Method (FEM). Therefore, SPH has drawn quite some attention in the geomechanics community [4,5,11,67]. Maeda et al. [38] and Naili et al. [48] are among the first to apply SPH to real geotechnical problems such as seepage analysis and soil liquefaction, but the constitutive models for soils they used were simple and not rigorous. Later, Bui et al. [5] incorporated the elastoplastic model based on the Critical State Soil Mechanics into the SPH framework to describe soil's behavior in large deformation and post-failure analysis. Since then, many researchers, following the SPH framework in [5], developed new formulations, numerical algorithms and implemented various constitutive models [2,49,53]. Plenty of reports can be reached in related areas including debris flow [7,14,52], landslide [8,13,51], slope stability [3,6,50], etc.
At the very beginning of those large deformation and post-failure problems is the emergence and development of shear strain localization, which is the key issue while modeling the granular material. However, it is well known that the classical continuum suffers from the pathological mesh-dependence [57]. To overcome this drawback, various regularization techniques have been developed by introducing an internal length into the constitutive equation [31], e.g. the nonlocal theory [60], strain gradient theory [64], micropolar theory [46] and theory with viscosity [1]. Among them the micropolar theory is chosen in this work because of its clear physical relevance for granular materials [18,26,59]. Micropolar theory, also known as Cosserat theory [20], considers the rotations as independent variables in addition to the translations in the classical continuum, which means that each material point has 6 DOFs (three for translation and the others for rotation). Besides, Micropolar theory also introduces some other physical quantities such as micro-curvature, defined as the gradient of the rotational angle, and couple stress energetically conjugate to the micro-curvature. As a result, the number of basic physical quantities is extended from 15 within the conventional continuum mechanics to 42 in the Micropolar theory (18 stress and strain components, 18 couple stress and micro-curvature components, and the aforementioned 6 DOFs). Therefore, it is foreseeable that the framework of the Micropolar theory can be of relative complexity in contrast to the conventional continuum. Nevertheless, the price is affordable as the advantages are distinct. For example, Micropolar theory is capable with estimating the width of the shear band due to the introduction of the internal length [59]. Moreover, this theory is well suited for granular materials such as sand, where particle rotations play an important role [63]. A comprehensive account of the theory can be found in the book by Vardoulakis [65].
As the extension of the classical continuum mechanics, micropolar theory has been applied to tackle various problems of progressive strain localization. Mühlhaus and and Vardoulakis [46] treated the shear band formulation as a bifurcation problem, and made theoretical prediction of the shear band thickness in agreement with experimental observation. Mühlhaus [47] analyzed the limit load of tunnel statics within the framework of micropolar theory by the finite element method. de Borst [15] presented the von-Mises elastoplastic model and analyzed the localization under static and dynamic loading [17], but this model is more suitable for metals rather than granular materials. For geomaterials such as stone and soil, de Borst [16] proposed the pressure dependent J 2 -flow theory for micropolar continuum, which is widely adopted in geomechanics [25,39,58]. Besides, many studies are conducted based on the micropolar hypoplasticity model [24,61,62]. Furthermore, the micropolar theory has been used in many other research areas, such as biomechanics [54], crystal plasticity [12,41] and composite material [9,55].
Despite the achievements mentioned above, it should be noted that the micropolar theory is no longer in force under pure tensile or compression as no rotational degrees exist in such circumstances. In addition, the micropolar theory is usually realized with the aid of standard finite element procedure, accounting for the limited applications to problems of small strain or weak discontinuity. This is where this paper steps in.
The motivation of this paper is to establish a general SPH framework based on the micropolar theory. As the first step, we consider an elasoplastic constitutive model based on the yield surface of Drucker-Prager. The performance of the Micropolar SPH (MPSPH) is compared with the conventional Euler kernel based SPH (CESPH) by conducting two typical geomechnics problems, i.e., biaxial compression test and sand column collapse. Meanwhile, a parameter analysis is given on the micropolar effect, defined as the influence induced by considering the microrotation. Simulation results show that the newly proposed method can deal with problems involving large deformation or discontinuous failure, and processes the capability in keeping the well-posedness of the boundary value problems with strain localization incorporated.
2 Micropolar continuum 2.1 Kinematics Figure 1 illustrates the three degrees of freedom in 2-D problems. Compared with conventional continuum mechanics which has only translational freedom, a point in the framework of micropolar continuum is also free in terms of rotation. This character leads to the asymmetrical strain rate tensor because the microscopic rotation is different from the macroscopic one, as shown in Eq. (1) In the equation, _ is the infinitesimal strain rate tensor, which is defined as where v denotes the velocity vector. w is the spin tensor with the expression w ¼ 1=2½rv À ðrvÞ T . Andis the microscopic rotation tensor, reading where f refers to the third-order permutation tensor and x ¼ _ h is the micro-rotation vector with h denoting the angle vector. Note that in a classical (Boltzmann) continuum, _ e coincides with _ as the slip and rotation between individual grains are neglected (i.e., -¼ w).
Apart from the strain tensor, there exists a curvature rate tensor in the micropolar theory, which describes the spatial derivative of the micro-rotation vector

Conservation laws
The conservation laws governing the physical phenomenon in micropolar continuum include the conservation of mass, linear momentum and angular momentum, i.e., and qJ dx dt À r Á l À qJc À f : where q is the mass density, r denotes the stress tensor, f the body force such as gravity in most cases, l the couple stress, c the body couple vector, and J the moment of inertia. Note that in this study we assume that the density is invariant, so Eq. (4) is not involved in simulation.

Cosserat elastoplastic model
The micropolar elastoplasticity is adopted to delineate the constitutive property of the Cosserat medium in this work, whose elastic part is shown as below in which _ r and _ l represent the stress rate and couple stress rate, respectively; _ e v ¼ trð_ eÞ with tr(Á) denoting the trace of a tensor; K; G andG c are the volumetric modulus, shear modulus and Cosserat shear modulus; and l denotes the characteristic length, which depends on the shape and the size of the micro-structure [17]. In fluid dynamics, the technique called Lagrangian velocities of tracer particles can be employed to measure the quantity experimentally [45]. In geomechnics, however, to the authors' knowledge, there is no reported method to determine the characteristic length. Theoretically speaking, the characteristic length of granular materials should be related with the averaged grain diameter, d 50 , since there is abundant experimental evidence that shear bands in granular media involve a significant number of grains [46].
To describe the plastic behavior for a Cosserat medium, the micropolar Drucker-Prager model is employed in this study, which has the following yield function where I 1 ¼ trr; c denotes the cohesion; k / and k c are the constitutive parameters related with the frictional angle, /, from Mohr-Coulomb model in which J 2 is the second deviatoric stress invariant in the micropolar continuum, which can be generalised as [16] J 2 ¼ a 1 s : s þ a 2 s : s T þ a 3 l : Herein s stands for the deviatoric stress defined as s ¼ r À 1 3 trðrÞI, where I is the identity tensor; a i is the material parameter with the requirement a 1 þ a 2 ¼ 0:5 to ensure that Eq. (11) can be retrieved to the conventional continuum form in the absence of couple stress. Notice that the set of values, a 1 ¼ a 2 ¼ 1 2 a 3 ¼ 1 4 , is used throughout the whole study [15,21].
By introducing the piecewise linear hardening/softening assumption into cohesion, we have, in which, c 0 is the initial cohesion; H denotes the hardening/softening paramter; and c represents the hardening parameter. In this paper, the micropolar equivalent plastic strain is chosen as the hardening parameter, owning the following rate form [16], with _ e p the plastic deviatoric strain. Besides, we adopt the non-associated flow rule with the following form of plastic potential where k w is another material parameter defined as with w the dilatancy angle.

Return-mapping algorithm
In this paper, we develop a one-step algorithm to determine the stress increment in a finite loading step, in which the gradient is evaluated at the trial stress state. It can be considered as the straightforward extension of the method proposed by de Borst et al. [15] for the pressure-independent Cosserat material. This algorithm begins with an elastic prediction, which leads to the trial stress, r t , and couple stress, l t , with the superscript t denoting the predicted value of a quantity. This practice is adopted throughout the whole work. Then the yield function Eq. (9) will be used to check whether the new stress state has reached the yield surface or not. If not, the updated stress and couple stress will be used as the input data for the next calculation loop. If the trial stress state goes beyond the yield surface, plastic failure occurs and a correction to the trial stress is needed. In this case, the corresponding plastic strain is computed with the plastic flow rule, i.e., Herein, r c and e p c represent the generalized Cosserat stress and plastic strain, with the former a block diagonal tensor incorporating both stress and couple stress while the latter assembling the plastic strain and curvature. Within the framework of Cosserat DP model, Dk, the increment of the plastic multiplier, has the following form, where c n is the hardening parameter at n-th calculation step. Eq. (17) features its simplicity but owns second order accuracy; details of the derivation are given in Appendix.
Unlike the classical FEM codes, where a global stiffness matrix is needed, the stress in the SPH framework can be updated directly once the increment of the plastic multiplier is determined [5]. One can acquire a straightforward cognition of the aforementioned whole process with the work flow shown in Fig. 2.

SPH formulations
SPH is a particle-based numerical method, in which physical information such as velocity and stress are carried by Lagrangian particles. The core of SPH is the kernel approximation and particle approximation.

SPH basics
A field function, f ðxÞ, and its spatial derivative, rf ðxÞ, at a given point x ¼ ðx 1 ; x 2 ; x 3 Þ can be approximated firstly with the following integral formulas: and rf ðxÞ h i¼ À where X is the supporting domain, W is the kernel function, and h is the smoothing length. In this study, the Wenland function [68] is used with the following expression in which, a d equals to 7=4ph 2 in 2-D and 21=16ph 3 in 3-D; q is non-dimensional distance between particles, expressed by q ¼ r=h, where r is the distance between two given particles.
With a further step, Eqs. (18)(19) can be rewritten by summing the contributions of all particles in the support domain of particle i, i.e., and where W ij ¼ Wðx À x 0 ; hÞ, n is the number of particles within the supporting domain of particle i, and m j =q j means the volume of the particle j. For more details of SPH, one can refer to Liu and Liu's book [32].

Discretization in MPSPH
By applying Eqs. (21)(22), Eqs. (1) and (3) can be discretized as where in which x ji ¼ x j À x i . Similarly, Eqs. (5-6) have the following forms: and Fig. 2 Flow chart of the return-mapping algorithm for Cosserat D-P model Compared with CESPH, except the widely-used Eq. (25), the other three formulas above are either additional or different in MPSPH discretization.

Artificial viscosity
Artificial viscosity, usually denoted by P ij , is a common numerical technique introduced to alleviate the unphysical oscillations in SPH method. Herein, the scheme proposed by Monaghan [43] is applied due to its simplicity and popularity: with in which, r ij is the vector pointing from j and i; a P and b P are two coefficients that need to be tuned according to specific problems; and c s is the artificial speed of sound, which will be further specified in the following section. In this study, the values of a P and b P are taken as 0.2 and 0, respectively, as such assignment can generate stable and reasonable results in all our simulations. Thus, incorporated with artificial viscosity, Eq. (25) is renovated as

Time integration
Among those available integration techniques [27], Verlet scheme [66] is chosen because of its low computational overhead, as it calculates the particle interaction only once per step. Thus, all the variables at each particle are updated with the following explicit integration: where n denotes the time step, and K traverses the following variables: velocity, v, angular velocity, x, stress, r, couple stress, l, strain, e, and curvature, j. In order to avoid the divergence of the integration through time, every N s time steps (N s ¼ 40 in this study), K should be calculated with a new expression Similar with other explicit integration algorithm, the time step is restricted with the so-called CFL condition, based on the acceleration term and the viscous diffusion term [23] Dt where v CFL is the CFL coefficient, a m ax is the maximum acceleration and / m ax, defined in Eq. (28), is the maximal term among all particles. In this work, the CFL coefficient is taken as 0.2 and c s ¼ 200 m/s is utilized. This will give rise to a sufficiently small time step, and thus, the computational stability can be achieved.

Boundary condition
Accurate enforcement of the boundary condition plays a crucial role in obtaining precise simulation results for any numerical method. Achieving this, however, is not an easy task for SPH as it suffers from the problem of particle deficiency near the boundary. Thus, various schemes have been developed to resolve the trouble including ghost particles [29], repulsive forces [44] and boundary particles [56].
In the present work, the solid boundary is modeled with three layers of dummy particles, which are either fixed or move with a prescribed velocity. To impose a non-slip condition, the boundary particles are reassigned with the following velocity and stress when interacting with soil particles: In these three equations, the subscripts b and s denote the boundary and soil particles, respectively. v w and a w , accordingly, are the wall velocity and acceleration. The superscripts a; b ¼ x; y; z refer to the Cartesian components. Note that Eq. (35) is only used for interaction, not the update of positions of boundary particles. Besides, the couple stress and angular velocity at all boundary particles are set equal to zero. This is a simplified way to deal with the rotational boundary conditions, but it is observed that reasonable results could be achieved, which will be shown in the subsequent section. In Sect. 4.2 where the simulations of biaxial compression are conducted, a method to handle stress boundary condition is required. The approach suggested by Zhao et al. [69] is adopted to model the flexible confined boundary condition, but a new criterion is proposed herein to determine which particles should be treated with this procedure. By doing so, Eq. (29) has a new expression: Dr 2h: where r c is the confining stress, and Dr represents the shortest distance between a given soil particle and the confined boundary. The new criterion implies that all those soil particles whose supporting domains are truncated by the confined boundary should take into account the effect of the confining pressure. Therefore, for an arbitrary soil particle, if Eq. (39) is satisfied, Eq. (38) will be applied, otherwise Eq. (29) is used in the calculation.

Numerical examples
In this section, the simulation of a simple shear test is conducted to validate the proposed SPH method in the first place. Subsequently, two typical geotechnical problems are considered to study the performance of MPSPH. First, the biaxial compression test is simulated to investigate the shear localization phenomenon, which is illustrated in detail with a parametric analysis on the micropolar effect (for definition, cf. Sect. 1). Then the proposed method is further examined with one of the most widely studied problems using SPH, i.e., the problem of sand column collapse. In both cases, results are compared with those from the conventional method, CESPH. In all simulations, the problem domain is discretized as particles in regular lattices with a particle spacing Dp. And the smoothing length h is fixed as ffiffi ffi 2 p Dp. All material parameters are summarized in Table 1. Note that Chen et al. [10] argue that J ¼ l 2 , but such a value will account for an unacceptable significant angular acceleration. Thus, we tuned the coefficient and set J ¼ 32 in all simulations except when discussing J-induced micropolar effect. Figure 3 illustrates the initial configuration of the squareshaped numerical model of 0.3 m length exploited in the simple shear test. The simulation model is divided into two separate parts, i.e., the centrally located specimen (100 mm Â 100 mm) and the boundary area enclosing the specimen. For realizing the pure shear condition, the velocity field is predefined with Eq. (40). Besides, the boundary particles are maintained with the initial state, while the specimen particles move freely during the whole simulation. The particle distance is taken as 10 mm, which leads to 961 particles in total. Three cases are conducted with values of confining stress of 50, 75, and 100 kPa.

Simple shear test
The stress path and stress-strain relations of the simple shear test with the micropolar Drucker-Prager model are plotted in Fig. 4. The simulations adopt the strain softening elastoplastic model, so the yield surface will shrink as the plastic strain develops. After shear starts, all loading stress paths move vertically upwards from their initial hydrostatic confined states, as shown by the hollow symbols in Fig. 4a. During the plastic flow, all stress states are enforced to remain on the yield surface. Thus the shear stress will be decreased in the softening regime, which is represented by the symbols with crossings in Fig. 4a and the declining lines in Fig. 4b. Besides, a good agreement is found between the numerical results and theoretical calculations, represented by the hollow symbols and the solid lines in Fig. 4b, respectively. A minor deviation is observed at the large strain stage. This is induced by the calculation error on the hydrostatic pressure, which can be verified with the post-failure stress states in Fig. 4a. In general, the results from the simple shear test indicate that MPSPH, incorporating the micropolar model and SPH, owns the potential to capture the appropriate stress state of the geomaterial with high accuracy. Figure 5 depicts the initial geometry and the boundary conditions in the plain-strain biaxial test. The soil sample is 50 mm wide and 100 mm high. Three resolutions, i.e., 2 mm, 3 mm and 4 mm, are chosen to study the shear localization problem. The characteristic length is taken as 1 mm. The techniques mentioned in Sect. 3.5 are applied to model the boundary conditions. In addition, the upper boundary particles move together downward at the speed of 10 cm/s, which can guarantee both accuracy and efficiency. The fixed and constant moving boundary conditions are well achieved by restricting particle displacement or velocity updates. One weak particle is introduced at the center of the sample to trigger the initialization of the shear localization. Besides, the specimen is confined with pressure 100 kPa.

Biaxial compression test
First, we simulate the biaxial test using CESPH with three numerical resolutions. Figure 6 portrays the equivalent strain contour at axial strain, e a , (ratio of the displacement of the top boundary to the sample height) equal to 6%, and equivalent strain is calculated as e eq ¼ ð 2 3 e Ã : e Ã Þ 1=2 , where e Ã ¼ À ðtr=3ÞI. As presented in Fig. 6, it is obviously observed that the X-shape shear band becomes wider with a larger particle spacing. In addition, a finer discretization can cause an increment in the maximum equivalent strain. These outcomes reflect the fact that the conventional continuum model suffers from the pathological mesh-dependence in localization analysis. By contrast, we find MPSPH exhibits less sensitivity to the change in numerical resolution. In Fig. 7a, all three resolutions lead to nearly the same shear band shape and size, and the close magnitude of the equivalent strain. Therefore, it is generally recognized that MPSPH gives rise to objective strain localization independent of numerical resolution. The size of shear bands is not a numerical result in MPSPH but is controlled by the characteristic length. This conclusion can be further validated with the rotation plot as shown in Fig. 7b. Hence, MPSPH clearly illustrates its capability of regularization in shear band pattern. Figure 8 presents the comparison of axial stress-strain relations between CESPH and MPSPH with different resolutions. Both methods predict that when axial strain reaches around 1%, the soil gets yielded and the peak stress is approximately 330 kPa. Starting from the same stress state, all curves experience almost identical elastic period before yielding. Following the yield stage is the strain softening plastic flow, which corresponds to the declining segments in the figure. In Fig. 8a, with finer resolution, the peak stress becomes slightly higher and more delayed. Additionally, the specimen displays more stiff behaviors in the softening regime of a coarse particle spacing than in a  fine discretization. In departure from the result within the framework of CESPH, the axial stress versus strain curves from MPSPH indicate that the mesh dependency problems have been remarkably alleviated. As illustrated in Fig. 8b, all curves with varied resolutions coincide with each other after peak values, until a small difference appears at the large deformation. Similar observations can be found in the existing literature [33,57]. Therefore, the regularization effectiveness of the proposed method, MPSPH, has been proven again, showing that not only the shear band patterns but the post-failure mechanical behavior is insensitive to the particle size.

Parametric analysis of micropolar effect
From Sect. 2, we can conclude that many coefficients play a part in the micropolar effect. In this section, we choose the following factors, i.e., the characteristic length, l, the Cosserat shear modulus G c and and the moment of inertia J, as they are typical parameters in micopolar theory. Firstly, the effect induced by the characteristic length is studied, as shown in Fig. 9. Note that all the data points are acquired from the dash line in the inset, which is a quarter width away from the right vertical boundary. The snapshot delineates the symmetric patterns for the equivalent strain   Fig. 7 Contours of a equivalent strain and b rotation for biaxial tests from MPSPH at e a ¼ 6% and the skew-symmetric patterns for the micro rotation, which is in line with expectation. Meanwhile, both the peak values of e eq and h decline as the characteristic length becomes larger, and the shear band gets wider. Similar patterns are observed in Fig. 10, which depicts the micropolar effect induced by the Cosserat shear modulus with a range from 0 to 4000 kPa. In case of G c equal to 0 kPa, the hollow square line portrays the result from CESPH without micropolar effect. In this condition, there is no microscopic rotation, so the corresponding line is not drawn in Fig. 10b. It is found that the peak values of both equivalent strain and rotation decrease following the increasing Cosserat shear modulus, which means micropolar effect gets stronger if the Cossearat shear modulus becomes larger. Meanwhile there exists a certain range beyond which the micropolar effect is not sensitive to the value change of Cosserat shear modulus. In this study, the range is between 1 and 1000 kPa.
The effect of the moment of inertia, J, is illustrated in Fig. 11. From the physical perspective, J describes the inertial resistance against rotational acceleration. This parameter can be ignored in quasi-static problems which is common in previous studies [19,28]. However, the moment of inertia is indispensable in MPSPH as it is an explicit method. Like the characteristic length, the moment of inertia is also dependent on the shape and the size of the micro-element. Herein, we explored J-induced micropolar effect by choosing several example values within a relatively large range. From Fig. 11, it is found that the moment of inertia is a key factor to influence the behavior of biaxial tests. The peak value of the equivalent strain drops as the moment of inertia increases. A similar tendency is observed for the microscopic rotation. Second, there exists a certain range (between 0.1 and 10 m 2 in this case), within which the biaxial test response is insensitive to the value change of the moment of inertia. Besides, in parameter analysis, an unrealistic result is observed when J is either too larger or small. Therefore, a calibration about the inertia of moment is recommended, before conducting simulations.

Sand column collapse
In this section, another numerical simulation, namely, sand column collapse, is investigated, which has been widely studied experimentally and numerically [11,35,36,40]. The initial geometry and boundary conditions are shown in Fig. 12. The column's initial width and height, d i and h i , are 200 mm and 100 mm, respectively. The sand column is discretized with two resolutions, i.e., 2 mm and 4 mm, so there are around 5000 and 1250 particles involved accordingly. The other configuration parameters is listed in Table 1. In order to highlight the micropolar effect, the artificial viscosity is switched off in this simulation.
To compare the performance of MPSPH and CESPH, the results from these two numerical methods as well as experimental data [35] are drawn in Fig. 13, which depicts the evolution of the runout distance. Those solid lines are obtained by recording the position of the column toe at each instant. Note that the dash lines represent the experimental process. Actually, in Lube et al. [35], there are only empirical formulas that describe the relationship between the incremental runout distance, dd, or collapsing duration, t 1 , and aspect ratio a ¼ h i =w i . Thus in dash lines, none but the turning points from curved lines to horizontal lines are plausible. However it is still meaningful to do the comparison between numerical and experimental results. All the simulations last 0.6 s and reach the stable stages finally. As shown in Fig. 13a, after about 0.5 s, the motion of the flow front ceases using CESPH and the incremental runout distance is around 0.274 m. However, in the case where l  Although these values from MPSPH vary with the data from the empirical expression from Lube at al. [35] (t 1 ¼ 0:33 s and dd ¼ 0:160 m), they are much closer to the experiment data than the result from CESPH. Moreover, when the characteristic length is assigned a smaller value, more precise result could be achieved. Similar conclusions could also be drawn from Fig. 13b where the mesh is coarser. Therefore, compared with CESPH, MPSPH can get more accurate results for different resolutions.
In addition, inside the sand column, there exists an interface which distinguishes the flow regime from the static domain. Traditionally, researchers usually use the strain contour to find such slip plane. Now, MPSPH could offer an alternative tool to search for the slip surface. Figure 14 demonstrates the contours for microscopic rotation and equivalent strain of the final deposit. The white dash lines are the so-called slip surfaces which separate the soils generally at rest from those after motion. From both images we can obtain right-angled trapezoid static regimes, which are similar with each other.

Closure
For the first time, the micropolar SPH is developed, which is based on the Cosserat Druck-Prager model considering strain softening. The one-step algorithm is proposed to update the stress state, with features in straightforward implementation and a second-order accuracy. Two numerical simulations are conducted to compare the performance of the newly proposed method with conventional SPH. Some advantages of the MPSPH are observed, including objective shear band patterns and stress-strain relations in the post-failure regime, independent of mesh size, and more accurate simulation without the employment of artificial viscosity. Besides, parametric studies are performed to investigate the influence of the micropolar effect. It is found that the characteristic length, Cosserat shear modulus and moment of inertia all play critical roles. This paper only reports our preliminary study on MPSPH. Several problems, i.e. better implementation of micropolar boundary condition, physical determination of micropolar parameters, require further investigation. Nevertheless, MPSPH shows to be a promising approach for geomechanical modeling as it can effectively regularize the strain localization in geomaterials.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.