Similar articles

Abstract: Standard finite element approaches are still ineffective in handling extreme material deformation, such as cases of large deformations and moving discontinuities due to severe mesh distortion. Among meshfree methods developed to overcome the ineffectiveness, Reproducing Kernel Particle Method (RKPM) has demonstrated its great suitability for structural analysis. This paper presents applications of RKPM in elasto-plastic problems after a review of meshfree methods and an introduction to RKPM. A slope stability problem in geotechnical engineering is analyzed as an illustrative case. The corresponding numerical simulations are carried out on an SGI Onyx3900 supercomputer. Comparison of the RKPM and the FEM under identical conditions showed that the RKPM is more suitable for problems where there exists extremely large strain such as in the case of slope sliding.

.  INTRODUCTION

 Geotechnical engineering’s well-developed finite element methods (FEMs) for geometrical and material nonlinearities facilitated much of the work accomplished. Nevertheless, standard finite element approaches are still ineffective in handling extreme material deformation, such as cases of large deformations and moving discontinuities due to severe mesh distortion. In order to overcome the ineffectiveness, some FEMs, such as Arbitrary Lagrangian Eulerian (ALE) method, which allow continuous remeshing during computation, were introduced (Hirt, 1974) and developed (Hughes et al., 1981; Belytschko et al., 1982; Donea, 1983). However, more effort is still required to go around the problem of the oscillation caused by the convective effect (Chen et al., 1996). To deal with these disadvantages, a new family of numerical methods was developed.

 All methods in this family, such as Smooth Particle Hydrodynamics (SPH) (Lucy, 1977; Monaghan, 1982), Particle in Cell Methods (PIC) (Harlow, 1964; Brackbill, 1986; Sulsky et al., 1994), Diffuse Element Methods (DEM) (Nayroles et al., 1992), Element Free Galerkin Methods (EFG) (Belytschko et al., 1994; Lu et al., 1994) and Reproducing Kernel Particle Methods (RKPM) (Liu et al., 1995; 1997), share a common feature in that no mesh is required and shape functions are formulated based on nodes, which distribute in the domain we want. More detailed classification and overview of meshfree methods, with their advantages and disadvantages, are presented in Fries and Matthies (2003).

 Among these methods, the EFG and the RKPM have been demonstrated as most suitable for structural analysis. Taking into account that the RKPM provides a general formulation for the construction of shape functions for meshfree methods, with specific discretization of the reproduced equation, so that the SPH and the EFG methods can be recovered. We present only the implementation of RKPM in this paper. Chen et al.(1997) extended RKPM to hyper-elasticity with large deformation.

 In the present paper, after a brief introduction to the RKPM, we present a meshfree discretization method for elasto-plastic problems with the Drucker-Prager and Mohr-Coulomb models that had been tested and proved suitable for geotechnical materials. A slope stability problem in geotechnical engineering is chosen as a sample case to demonstrate the advantages of the meshfree methods.

.  REVIEW OF THE RKPM

 We begin with the integral transformation of function u(x) by (1)

 in which u^a(x) is called as the reproduced function of u(x), and the modified kernel function, that is expressed as (2) (3) where a is the dilation parameter of kernel function Φ_a(x−s), and C(x; x−s) is called as the correction function, which is constructed to avoid the difficulties resulting from finite domain effects and to minimize the amplitude and phase error. The correction function is expressed by an Nth order polynomial in (x−s), that is (4) where (5) (6)

 and b_i(x) are functions of x, which are determined through satisfying the reproducing conditions. From Taylor series expansion, one can obtain (7) where Finally, the modified kernel function is obtained as (8) where The discretized reproducing equation is obtained by performing numerical integration in Eq.(1). By discretizing the reproducing equation with the trapezoidal rule, we get: (9)

 and (10) where n is the total number of particles, and can be interpreted as the shape functions of u^a(x). In principle, we take Φ_a(x−x_i) to be positive and set its maximum value at x=x_i. The function should rapidly approach zero after |x=x_i| exceeds a small value so that the shape function associated with node i has an influence on only a small group of surrounding nodes to provide computational efficiency.

 The extension of Φ_a to a multi-dimensional case can be achieved by using tensor products of the one-dimensional kernel functions (11)

 In this construction, the supports of kernel functions are rectangular and hexahedral in geometries, each with a center x_i and dimensions 2a₁×2a₂ and 2a₁×2a₂×2a₃, in two- and three-dimensional problems, respectively.

 The correction function for the multi-dimensional case is expressed as (12) where The coefficients are determined from the reproducing conditions. Finally, the multi-dimensional RKPM interpolation is obtained as (13) where (14)

 and ΔV_i is the volume associated with particle i.

.  RKPM FOR ELASTO-PLASTIC PROBLEMS

 Different from those in the FEM, the shape functions in the RKPM do not possess the Kronecker delta properties, i.e., . This will cause difficulties in numerical implementation while dealing with kinetic constraints on the geometrical boundaries of a model. A modified RKPM shape function that possesses Kronecker delta properties is required for numerical implementation.

 Recall that (15) (16) where By substituting Eq.(16) into Eq.(15), one can obtain (17) where Note that (18)

 and (19)

 are nodal values. u^a and δu^a satisfy the following boundary conditions: (20) where denotes a particle numbers set in which the associated particles are located on From Eq.(18), we directly obtain: (21)

 By using the expressions obtained above, we can discretize the weak forms of the elasto-plastic problems. As we known well, the equilibrium equation of the small strain elasto-plastic problem can be written as (22)

 Its weak form based on virtual displacement principle can be written as (23)

 with the constitutive relations or (24)

 Similar to those in finite element methods, the following expressions must hold theoretically (25)

 in which B^T and are defined as u= δu= ε= δε= respectively. B has a well known relation with as (26)

 in which L is a differential operator between strains and displacements. However, different from that in finite element methods, volume integration of Eq.(25) over the domain Ω is not the sum of the individual elements, but the points.

 In addition, though Eq.(25) hold theoretically, it will not generally be satisfied at any stage of numerical computation, owing to the nonlinearities of elasto-plastic problems. Namely, the residual force vector is (27)

 for evaluating the material tangent stiffness matrix.

 At any stage, the incremental form of Eq.(25) has to be employed. This implies that, for an increment of load, we have (28)

 Based on this equation, the numerical procedure, such as Newton-Raphson Method and Tangential Stiffness Method (sometimes named as Generalized Newton-Raphson Method), can be directly applied.

.  CASE STUDY

 For the frictional materials such as soil and rock in geotechnical engineering, a non-associated or associated Drucker-Prager (and Mohr-Coulomb) plasticity model is appropriate for modeling its constitutive behavior. This can be described generally as follows (Drucker and Prager, 1952).

 A quadratic stored energy function Ψ(ε^e, ξ), which results in linear elasticity and linear hardening, is defined as (29) where c^e and H are constant modulus tensors. Note that Ψ(ε^e, ξ) could be defined such that nonlinear elasticity and/or hardening would have resulted. The fourth-order tensor c^e is a tensor of isotropic elastic tangent modulus defined as (30) where is the elastic bulk modulus, and are the Lame parameters, (1)_ij=δ_ij is the Kroneker delta, and (I)_ijkl=(δ_ikδ_ji+δ_ilδ_jk)/2 is the fourth-order identity tensor.

 Let the strain-like vector of plastic internal variables ξ : Ω×[0,T]→(² have a volumetric component and a deviatoric component (31) where and The hardening/softening modulus matrix H is defined as (32) where and are the bulk and shear hardening/softening modulo, respectively. The stress-like vector of plastic internal variables α is defined as (33)

 A Drucker-Prager yield function takes the form (34)

 with derivatives where and , in which Note that for the continuum mechanics convention used throughout this paper, compression implies p<0. The above material constants and β may be defined in terms of the cohesion and friction angle used to describe a Mohr-Coulomb material (Owen and Hinton, 1980) (35)

 The value A=−1 coincides with a cone that circumscribes the Mohr-Coulomb envelope−passing through its outer apexes−in three-dimensional stress space, and A=1 coincides with a cone that passes through the inner apexes of the Mohr-Coulomb envelope. Because of the form of in Eq.(34), standard bulk and shear hardening/softening through α₁ and α₂ cause the size of the yield cone to change, and not its shape. Incorporating additional hardening/softening through β would allow frictional hardening/softening that is appropriate for modeling a cohesionless granular material, such as sand, thus the change of slope of the yield cone would be possible.

 A plastic potential function φ(σ, α) is defined similar to the yield function ϕ(σ, α) as (36)

 with derivatives (37) where b is the material dilation constant. Associated plasticity results if β=b, but typically for soil and rock, this is not the case. Usually there exists β>b with b>0 for a dilatant material, and b<0 for a contractant material. By setting β=b=0, the J2 flow (von Mises) plasticity model is recovered, which is useful for modeling the undrained conditions in a cohesive soil.

 With the plastic potential function φ(σ, α) as defined in Eq.(36), the evolution of ξ then becomes (38)

 Note that b takes the form (39)

 which is analogous to the dilatancy factor used by Rudnicki (1977) and Rudnicki and Rice (1975).

 There are many kernel functions proposed. Here we employ a cubic spline function as the kernel function of Eq.(11) (40) where (41)

 As an illustrative case, we study a plane strain problem of slope stability in geotechnical engineering. The geometrical model of the slope is shown as Fig.1 with a slope of 5/3 and is discretized with 3029 points (Fig.3). The constitutive relation is Mohr-Coulomb model. Under self-gravity, the slide will happen while the friction and cohesive strength of soil decreases.







 During the numerical computation, the cubic spline function is chosen as kernel function Eq.(40), with the cohesive-like strength parameter of 0.015, friction-like parameter of 0.04, dilation and hardening parameter of zero. As a result, the maximum relative and maximum absolute errors are 10⁻¹⁰ and 10⁻¹², respectively. The deformed shape of the slope is depicted in Figs.4 and 6. In order to compare with mesh-based methods, such as the finite element method, the computed results by the meshfree method is shown with and without a mesh as Figs.5 and 2, respectively. The extremely deformed part near the slope-toe is outlined with a mesh in Figs.7 and 8.

.  CONCLUSION

 The numerical computation shows that, while the deformation is not sufficiently large, such as in an elastic stage, the FEM and the RKPM will produce almost the same results, but the meshfree method costs more computation time. However, while the deformation is large enough, even with the same discretization and allowable errors, the FEM would result in divergence with RKPM resulting in convergence. As to the situation of extreme deformation, for example the case of slope sliding we are discussing here, the FEM cannot deal with it due to the non-positiveness of the Jacobian, which is caused by the severe mesh distortion. The conclusion can be drawn that for extreme deformation, the meshfree method RKPM is more appropriate.

 In addition, it should be mentioned that parallel computation was used for complete comparison between FEM and RKPM on an SGI Onyx3900 supercomputer with 1, 2, 4, 8 and 16 processors utilized.