On the application of the maximum entropy meshfree method for elastoplastic geotechnical analysis

Computers and Geotechnics 84 (2017) 68–77

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

On the application of the maximum entropy meshfree method for elastoplastic geotechnical analysis Omid Kardani a, Majidreza Nazem b,⇑, Mina Kardani c, Scott Sloan d a

School of Engineering, Deakin University, Melbourne, Australia School of Engineering, RMIT University, Melbourne, Australia c Software Engineer, ITW Construction, Asia Pacific, Melbourne, Australia d Australian Research Council Centre of Excellence for Geotechnical Science and Engineering, The University of Newcastle, Australia b

a r t i c l e

i n f o

Article history: Received 18 May 2015 Received in revised form 12 October 2016 Accepted 22 November 2016

Keywords: Maximum entropy meshfree method Dynamic relaxation Analytical stress integration Elastoplastic geotechnical analysis

a b s t r a c t In this study, the Maximum Entropy Meshfree (MEM) method is employed for analysing geotechnical problems involving material nonlinearity, assuming small strains. The efficiency of the MEM method is evaluated through several solution schemes for the global governing equations as well as the local constitutive equations. The conventional implicit approach involving the Newton-Raphson method and an explicit adaptive dynamic relaxation technique are employed for solving the governing equations, while local constitutive equations are solved numerically as well as analytically. Two- and three-dimensional numerical experiments are performed to study the efficiency of different configurations of the solution scheme, which leads to some important conclusions about application of the MEM method in geotechnical problems. Crown Copyright Ó 2016 Published by Elsevier Ltd. All rights reserved.

1. Introduction Meshfree methods have undergone significant progress during the last two decades. The key idea behind these methods is to utilise a set of arbitrarily-distributed particles without any mesh to represent the state of a system and to record its movement. Therefore, meshfree methods have the potential to be superior to conventional mesh-based methods, such as the finite element (FE) method, in the sense that they do not suffer from computational difficulties resulting from mesh structure such as mesh distortion, complex mesh generation, and remeshing. This feature has made meshfree methods an attractive choice for geotechnical engineering problems involving large deformation, discontinuities, and moving boundary conditions. The Element-Free Galerkin (EFG) method is a widely-used meshfree method in computational geomechanics, thanks to its robustness, relatively high accuracy and superior convergence [22,28,32,34]. The EFG method in its original format as developed by Belytschko et al. [6], exploits Moving Least Square (MLS) shape functions to interpolate the field variables. However, as the MLS functions do not satisfy the Kronecker delta property, their application may give rise to problems in imposing the essential boundary

⇑ Corresponding author. E-mail address: [email protected] (M. Nazem). http://dx.doi.org/10.1016/j.compgeo.2016.11.015 0266-352X/Crown Copyright Ó 2016 Published by Elsevier Ltd. All rights reserved.

conditions. As a remedy, alternative shape functions and, consequently, several meshfree methods have been developed for analysing problems in solid mechanics and geomechanics. Examples include, but are not limited to, the Point Interpolation Method [17], the Radial Point Interpolation Method [35,13], and the Local Radial Point Interpolation Method [18]. The major difference among these methods is the type of shape functions employed for interpolating a field variable. Recently, the concept of maximum entropy (max-ent) shape functions, originally developed by Sukumar to formulate polygonal interpolants [30], was introduced to the EFG method [3,31]. Due to its robustness and computational stability, the EFG method with max-ent shape functions, or otherwise called the Maximum Entropy Meshfree (MEM) method, has been a popular numerical framework for analysis of engineering problems, including incompressible elasticity [23], nonlinear analysis of reinforced concrete structures [26], analysis of thin shells [19] and coupled FE-MEM analysis of linear and nonlinear problems of solid mechanics [32]. Nevertheless, the application of MEM to geotechnical engineering is relatively new. Recently, Nazem et al. [20] employed this method for consolidation analysis of porous media. In this study, the application of MEM to incremental elastoplastic analysis of geotechnical problems is investigated. To present a comprehensive evaluation, several approaches are utilised for solution of the governing and local constitutive equations.

69

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

In this study the governing equations are solved in an incremental manner by the Newton-Raphson iterative solution scheme or, alternatively, by an explicit time integration approach. The former is a widely used Newton-Raphson method for solution of nonlinear system of equations and benefits from unconditional stability as well as efficient convergence properties. However, being an implicit method, it requires the formation and solution of a linear system of equations in each iteration. An alternative approach is to undertake an explicit approach to avoid dealing with linear systems of equation. This is of particular interest where large scale problems are addressed, as in 3D analysis of geotechnical problems. For this purpose, adaptive dynamic relaxation is proposed in this study. Furthermore, the elastic-perfectly plastic von Mises model is considered for the constitutive relations. The solution of local constitutive equations is also performed by two different approaches, i.e. numerically and analytically. The numerical method used to deal with stress-strain relations is the backward Euler method, while analytical solutions are obtained using a return mapping scheme. The structure of this paper is as follows. In Section 2, the max-ent approximation is briefly discussed and the governing equations for small strain incremental elastoplastic analysis within the framework of the MEM are derived. Then, implicit and explicit approaches for solution of the governing equations are presented in Section 3, followed by a study of numerical and analytical stress integration schemes. Numerical results on two- and three-dimensional footing problems are presented and discussed in Section 5. Finally, in Section 6, key findings of this study are summarised and conclusions are drawn.

The MLS shape functions do not satisfy the Kronecker delta property. This can give rise to difficulties in applying the essential boundary conditions. With this in mind, a class of interpolating functions that satisfy the Kronecker delta property, known as the Maximum Entropy (max-ent) shape functions, was recently introduced to the EFG framework to facilitate imposition of essential boundary conditions [3,31]. 2.1. Max-ent shape functions The principle of maximum entropy was first introduced by Jaynes [11]. According to this principle, the most likely distribution of discrete probabilities is obtained by solving n X Maximise : HðPÞ ¼  Pi ln Pi i¼1 n X subject to : Pi ¼ 1; i¼1

Background cell

i¼1

where P i  Pðvi Þ; i ¼ 1; . . . ; n are the probabilities of the occurrence of n discrete events vi, H represents the entropy of a discrete probability distribution, and Eðg r ðvÞÞ is the mathematical expectation of a function g r ðvÞ. Based on the analogy between the discrete probabilities in (1) and shape function values, max-ent shape functions were first proposed as polygonal interpolants [30]. The local max-ent functions were then introduced in [3], which are more attractive for meshfree methods as they can be defined for a local set of nodes. Weight functions were later introduced as a compact support to max-ent shape functions [3,31]. In the presence of weight functions, wi, the max-ent shape functions, Ui, correspond to the solution of the following constrained optimisation problem

2. Element-Free Galerkian

HðU; wÞ ¼ 

Maximise :

The Element Free Galerkin (EFG) method is a meshfree method which is based on a weak form of the global governing relations [6]. In its original form, the EFG method utilises Moving Least Square (MLS) shape functions to approximate the field variables and the essential boundary conditions are imposed using Lagrange multipliers. Moreover, a background mesh is employed only for the purpose of numerical integration. In this way, a two dimensional continuum with problem domain X bounded by domain boundary C is discretised by the EFG method as shown in Fig. 1. Similar discretisations can be applied for three-dimensional problems.

subject to :



n X

Ui ln

n X

i¼1 n X

i¼1

i¼1

Ui ¼ 1;

Ui



wi

Ui xi ¼ x;

n X

i¼1

i¼1

Ui yi ¼ y;

Ui zi ¼ z

and are derived as

Ui ¼

wi ek1 ðxi xÞk2 ðyi yÞk3 ðzi zÞ

n X

wj e

ð3Þ

k1 ðxj xÞk2 ðyj yÞk3 ðzj zÞ

j¼1

y

Support domain of an integration point

Domain boundary Γ

z Traction boundary ΓT

Essential displacement boundary Γu Field node

Domain boundary Background mesh

n X

ð2Þ

Surface traction

Integration point

ð1Þ

n X P i g r ðvi Þ ¼ Eðg r ðvÞÞ

Problem domain Ω

Influence domain of a Field node

Fig. 1. Problem discretisation by the EFG method.

x

70

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

where x, y, and z represent the spatial coordinates of a field node, k1, k2, and k3 are the Lagrange multipliers associated with constraints in Eq. (2) and are obtained by solving the following convex minimisation problem using Newton’s method [31]

ðk1 ; k2 ; k3 Þ ¼ argmin ln

n X wj ek1 ðxj xÞk2 ðyj yÞk3 ðzj zÞ

!

ð4Þ

Upon discretising the problem domain with a set of nodes, an approximation of the displacements at point x, i.e. uh ðxÞ, can be obtained using the max-ent shape functions for a field node i at a point x, represented by Ui , as follows

uh ðxÞ ¼ ½ ux

uy

uz T ¼

j¼1

if

ri 6

if

1 2

if

1 2

< ri 6 1

ð5Þ

ri > 1

mi

from the ith field node to the point of interest and dmi ¼ as ci is the influence domain of the ith node defined, where as is a scalar parameter and ci represents the distance between two adjacent nodes in a uniformly distributed set of nodes. For a non-uniform nodal discretisation, ci is obtained by calculating an average nodal spacing in the influence domain of the ith field node. Effective values of as are between 2.0 and 4.0 [16]. Based on numerical results presented by Nazem et al. for geotechnical problems [20], a typical value of 3.0 is adopted in this study. Circular and spherical support domains are adapted for two- and three-dimensional problems in this study, respectively. 2.2. Governing equations In this section the governing equations for three-dimensional static analysis of geotechnical problems are presented. Note that the modification of the following equations for a twodimensional analysis is straightforward. Consider a threedimensional continuum with the problem domain X bounded by C. Equilibrium is satisfied provided that

rT r þ b ¼ 0 in X

ð6Þ

where r ¼ ½ rxx ryy rzz rxy ryz rxz  is the Cauchy stress vector, b denotes the body force vector, and r is the following differential operator T

@ @x

0

0

@ @y

0

r¼6 60

@ @y

0

@ @x

@ @z

0

0

@ @z

0

@ @y

6 6 4

@ @z

3T

7 7 07 7 : 5

ð7Þ

The boundary conditions associated with equation (6) are given

Z

X

Z

X

Z

duT b dX 

BT Dep BdX Z

Fint ðuÞ ¼ Z Fext ¼

X

X

ð12Þ

BT rðuÞdX;

Ui b dX þ

Z CT

Ui T dC

ð13Þ

ð14Þ

in which B is the matrix of derivatives of shape functions and Dep is the material constitutive matrix. A key issue in meshless methods is how to evaluate the integrals in Eqs. (12)–(14). The integration strategies for computing these matrices and vectors in the global equations are not unique, and are basically divided into two main categories; cell-based integration and nodal integration. The former, which uses the field nodes as integration points, does not require any cells (background mesh) and, consequently, provides a truly meshless approach. However, due to its spatial instability, the nodal integration strategy needs an extra stabilisation procedure [5]. In addition, the applicability of this strategy within the Maximum-Entropy framework has rarely been addressed in the literature. The cell-based integration approach, on the other hand, is more popular and has been widely used for integrating the global equations in the Maximum-Entropy Meshless method. To avoid further complexity, this approach is employed here for integration. It is notable that, for geotechnical applications, this strategy provides a stable solution [20,32]. 3. Pseudo-time integration In an incremental step-by-step solution of Eq. (11), it is assumed that the nodal displacements t u and internal nodal forces F at time t (load increment t) are known. Then, considering an appropriate time increment Dt, the nodal displacements at time t þ Dt, tþDt u, are calculated by solving

F

ð8Þ

 is the vector of prescribed where u is the displacement vector, u displacements on the displacement boundary Cu , n denotes the outward unit vector normal to the boundary C and T represents the prescribed traction on the traction boundary CT . The weak form of equations (6) may be written as

dðruÞT r dX 

ð11Þ

where Kep is the elastoplastic stiffness matrix and Fint and Fext denote the internal and external force vectors, respectively. The components of (11) are given by the following expressions

tþDt ext

Essential Boundary Conditions :rn ¼ T on CT  on Cu Natural Boundary Conditions :u ¼ u

X

Fint ðuÞ ¼ Kep u ¼ Fext

t int

@ @x

as

Z

ð10Þ

in which n is the number of field nodes in the support domain of point x, ui are nodal values of displacement, and ux , uy and uz are the displacements in the x, y and z directions, respectively. After simplification, the equilibrium equation is expressed by

Kep ¼

in which r i ¼ ddi is the normalised distance, di denotes the distance

2

Ui ui

i¼1

Various weight functions may be employed in the formulation of max-ent shape functions [26] as long as they are positive within the support domain, zero outside the support domain, and monotonically decreasing, yet sufficiently smooth near the boundary of a support domain. In this study, the following weight functions based on 2nd order continuous cubic spline approximation are used due to their popularity [34,36]

8 2  4r 2i þ 4r3i > > >3 < wðr i Þ ¼ 43  4r i þ 4r2i  43 r 3i > > > : 0

n X

CT

duT T dC ¼ 0:

ð9Þ

 tþDt Fint ¼ 0

ð15Þ

where the left superscript refers to the time increment. In an elastoplastic analysis, Eq. (15) is actually a nonlinear system of equations since the internal forces depend nonlinearly on the nodal displacements. This system can be solved in an implicit or explicit manner. 3.1. Implicit approach An efficient method to solve the nonlinear systems of Eq. (15), which is also widely used in nonlinear finite element analysis, is the Newton-Raphson scheme. This utilises an iterative approach

71

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

to solve system (15) based on the following set of relations for n ¼ 0; 1; 2; . . . tþDt

tþDt ext Kep F  tþDt Fint ðnÞ ðnÞ Duðnþ1Þ ¼

ð16Þ

tþDt

uðnþ1Þ ¼ tþDt uðnÞ þ Duðnþ1Þ

ð17Þ

t int with the initial conditions tþDt uð0Þ ¼ t u and tþDt Fint ð0Þ ¼ F . The termi-

nation criterion for the iterative scheme in (16) and (17) is based on the ratio of the norm of the unbalanced forces as

ktþDt Fext  tþDt Fint k2 6e ktþDt Fext k2

ð18Þ 3

6

in which the tolerance e is typically chosen in the range 10 to 10 [16]. Note that relation (16) is, in fact, a linear system of equations. This means that, in each iteration of the Newton-Raphson method, the global stiffness matrix must be formed and factorised to solve (16). Although iterative solution schemes may be employed to solve (16), they are not usually as effective as a direct solution based on Gaussian elimination unless very large and relatively well-conditioned systems of equations are involved [4].

ð19Þ

where the superscript dots denote the order of time derivatives, and M and C are fictitious mass and damping matrices, respectively, which are chosen such that convergence rate of the dynamic problem toward the steady-state solution is accelerated. Employing a diagonal mass matrix produced by the mass lumping technique and a diagonal Rayleigh mass proportional damping has some significant computational advantages. The assembly and factorization of global matrices are avoided due to uncoupled system of algebraic equations in which the solution component maybe computed independently. In addition, Belytschko and Mullen [7] reported that the use of the lumped mass technique leads to lower estimates of the linear system eigenvalues and can, therefore, increase the numerical stability limit of the explicit time integrator. Above all, it has been noted by Krieg and Key [14] that the accuracy of the explicit time integrator is improved as a result of the errors introduced by lumped masses and the central difference operator being compensatory. The Rayleigh mass-proportional damping adopted in this study is of the following form [21]:

c2R

v ðnþ1=2Þ ¼

ð20Þ

Other choices for mass and damping matrices have also been studied in various applications [2,27].After constructing the fictitious dynamic problem, its solution is obtained through an explicit

1 ðuðnþ1Þ  uðnÞ Þ h

ð21Þ

1 ðv ðnþ1=2Þ  v ðn1=2Þ Þ; h 1 ¼ ðv ðnþ1=2Þ þ v ðn1=2Þ Þ 2

v_ ðnÞ ¼ v ðnÞ

where h is a fixed time step. By substituting Eqs. (20) into (21), the governing equations for advancing the velocity and displacement can be obtained from the following set of relations for n ¼ 0; 1; 2; . . .

tþDt

The formation and solution of linear system of equations is the computational bottleneck when dealing with large scale threedimensional problems. In order to avoid this, the time integration can be performed in an explicit manner. An efficient method for this purpose, which makes it possible to use explicit dynamic algorithms to solve a static problem, is called dynamic relaxation (DR) [33]. DR is based on two basic observations. Firstly, the solution to a static problem can be seen as the steady-state solution of a dynamic problem. Secondly, when we seek the steady-state solution, the dynamic components, including mass, damping and the time step, are of no interest [33]. In this way, instead of solving the static problem (15) for each load increment, we formulate a fictitious dynamic problem whose steady-state response matches that of the original static problem. The fictitious dynamic problem for a desired load increment is given by:

C ¼ cM;

ext Mv_ ðnÞ þ CvðnÞ þ Fint ðnÞ ¼ F

v ðnþ1=2Þ ¼

3.2. Explicit approach

€ þ Cu_ þ Fint ¼ Fext Mu

time integration method. The central difference technique is an efficient choice as it is very simple to implement and also benefits from the largest stability limit among the second-order accurate integration schemes [15,21]. Besides, Oakley and Knight [21] proposed the use of a half-station central difference method as it possesses a smaller order of truncation error OðDt 2 =4Þ compared with OðDt 2 Þ for the full-station scheme. Thus, after the time discretisation of problem (19) and introducing the velocity vector v , the fictitious dynamic equations for the nth time increment can be written as:

    2  ch 2h M1 ðFext  Fint v ðn1=2Þ þ ðnÞ Þ 2 þ ch 2 þ ch

uðnþ1Þ ¼ uðnÞ þ hv nþ1=2

ð22Þ

Note that applying the matrix M1 is trivial as M is a diagonal matrix. Moreover, the velocities for the first time step can be calculated from

vð1=2Þ ¼

    2  ch h M1 ðFext  Fint ðuð0Þ ÞÞ v ð0Þ þ 2 2

ð23Þ

with v ð0Þ and uð0Þ being the initial velocities and displacements, respectively. The same convergence criteria as in (18) can be employed for the DR scheme as well. To achieve the optimum convergence rate for (22), the fictitious mass matrix M can be approximated based on the CourantFriedrichs-Lewy conditions for stability of the central difference scheme [8] as

Mii P

2 m h X jKep j 4 j¼1 ij

ð24Þ th

where Mii is the i diagonal entry of M, m is the number of columns in the global elastoplastic stiffness matrix Kep , and h is an arbitrary fixed step size (usually taken to be the unit step size) [21]. In addition, the fictitious damping coefficient c can be estimated based on mass-stiffness Rayleigh quotient as

v Tðn1=2Þ DðnÞ v ðn1=2Þ vTðn1=2Þ Mvðn1=2Þ

c¼2

!1=2 ð25Þ

in which DðnÞ is the diagonal estimator of the directional stiffness after the nth step and is given by





DðnÞ

 ii

   Fint  Fint ðnÞ ðnÞ i   i ¼ h v ðn1=2Þ i

ð26Þ th

where the outer subscript i denotes the i vector element. Finally, for highly nonlinear problems, the stiffness may undergo significant changes during the course of analysis. In order to accommodate these changes and ensure stable convergence of the method, the integration parameters need to be updated in an adaptive manner. The damping coefficient can be computed cheaply in each DR iteration. However, the calculation of the mass

72

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

matrix through (22) is usually computationally expensive due to the required stiffness calculation. It is, therefore, desirable to update the mass matrix only when the stability of the algorithm is critically perturbed. Based on this idea, Park and Underwood [24] proposed an adaptive dynamic relaxation (ADR) scheme where the fictitious mass matrix M is updated using (24) whenever the ADR algorithm becomes potentially unstable. In this regard, they proposed the following stability criteria based on the perturbed apparent frequency error estimators ei which is given by

ei ¼

2 h jðv_ ðnÞ Þi  ðv_ ðn1Þ Þi j : 4 jðuðnÞ Þi  ðuðn1Þ Þi j

ð27Þ

The ADR algorithm based on the above guidelines has been employed in finite element analysis of some geotechnical problems [10,12]. However, the application of ADR in the framework of meshfree methods has not been studied for problems in geomechanics. 4. Constitutive relations When the deformations are assumed to be infinitesimal, the strain-displacement relation is given by

e ¼ Bu

ð28Þ

in which u and e are the displacements and strains vectors, respectively, and B is the strain-displacement matrix. In an elastoplastic analysis, the stress increment corresponding to a total strain increment needs to be determined by integrating the following differential equation [36]

r_ ¼ Dep ðrÞe_

ð29Þ

in which superior dots indicate derivatives with respect to time and Dep ðrÞ is the elastoplastic constitutive matrix defined as

 T D @F @F D Dep ðrÞ ¼ D  @rT@ r : @F D @@Fr @r

ð30Þ

where D denotes the elastic constitutive matrix and F is the yield function. In this study, we employ the elastic-perfectly plastic von Mises model to predict the soil behaviour under an undrained condition. The yield surface in the von Mises model is defined by

FðrÞ ¼

1 T 2 s sk 2

ð31Þ

in which k is the yield stress in pure shear and s represents the deviatoric stresses given by

s ¼ r  mp with

p

ð32Þ

being

the

mean

stresses,

i.e.

p¼

1 mT 3

r,

and

m ¼ f1; 1; 1; 0; 0; 0gT . Therefore, with the displacements uðnÞ known, the stresses in the nth time increment are calculated at each integration point by

@F @r k P 0;

BðuðnÞ  u0 Þ ¼ D1 ðrðnÞ  r0 Þ þ k FðrðnÞ Þ 6 0;

kFðrðnÞ Þ ¼ 0;

ð33Þ

where k is a plastic multiplier and the subscript 0 refers to the initial state at the beginning of the physical load step. The stress rðnÞ is int then utilised in (13) to calculate Fint ðnÞ ¼ F ðuðnÞ Þ, which is then used

in (16) or (22) to obtain uðnþ1Þ . 4.1. Numerical approach For a given strain increment De ¼ BDu, the constitutive relations (29) must be integrated at each Gauss point. A pseudo time

stepping scheme can be used to integrate Eq. (29) numerically. In this way, the stresses are calculated at time t0 þ Dt based on the known stresses at the starting time t 0 . Moreover, in the static displacement analysis, the strain rate is assumed to be constant and given by

e_ ¼

De Dt

ð34Þ

0 By defining the pseudo time parameter T ¼ tt where Dt 0 6 T 6 1, and making use of (34), (30) can be differentiated as

dr @F : ¼ Dep De ¼ DDe  DkD @r dT

ð35Þ

Eq. (35) is a differential equation with initial value, which is the known stresses at time T ¼ 0, i.e. r0 . The numerical method used in this study to solve the initial value problem (35) is based on the stress integration scheme proposed by Sloan et al. [29]. It employs a combination of explicit Euler and modified Euler procedures, together with an automatic sub-incrementing scheme applied to the imposed strain increment De, to control the error in approximate integration of the constitutive relations (see [1] for details). 4.2. Analytical approach Given a strain increment De ¼ BDu, the local relations (29) can be solved analytically through a return mapping scheme. In this way, first the trial elastic stresses, rtr , are calculated as

rtr ¼ DDe

ð36Þ

Next, depending on the state of trial stresses, one of the following cases may happen which is dealt with accordingly to calculate the stresses in the time increment – if Fðrtr Þ 6 0, the material is in the elastic state and rðnÞ ¼ rtr . – if Fðrtr Þ > 0, the material has shown plastic response, so the new stress state is calculated as

rðnÞ

sffiffiffiffiffiffiffiffiffiffi 2 2k ¼ ptr þ str T str str

ð37Þ

where ptr ¼ 13 mT rtr and str ¼ rtr  mptr . 5. Numerical examples The performance of the Maximum Entropy Meshless method with different solution schemes is investigated by considering the classical geotechnical problem of the bearing capacity of an undrained layer of soil under a rigid rough footing. Both twoand three-dimensional analyses are performed. Fig. 2 shows the problem domain, boundary conditions and soil properties for the two- and three-dimensional analyses. The soil layer is assumed to be weightless and its nonlinear behaviour is predicted by the von Mises material model, assuming small deformations. Two different grids are used for analysing the footing problems, including a relatively fine and a relatively coarse grid. The properties of the grids including the number of field nodes, total number of Gauss points and the number of background cells are summarised in Table 1. Finally, the configuration of the solution scheme varies by employing explicit and implicit approaches for time and stress integration. In this regard, the following notation is utilised for various underlying procedures in the solution schemes: NR: refers to implicit time integration where the NewtonRaphson method is used. ADR: refers to explicit time integration where Adaptive Dynamic Relaxation is used.

73

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

B

Smooth boundary

Undrained Young’s modulus, Eu = 100 kPa Poisson’s ratio, vu = 0.495 Undrained shear strength, su = 1 kPa Undrained friction angle,

φ

u

5B

=0

B=1m

Smooth boundary

B

10 B

y

x Smooth boundary

(a) Two-dimensional analysis B

B Undrained Young’s modulus, Eu = 100 kPa Poisson’s ratio, vu = 0.495

5B

Undrained shear strength, su = 1 kPa Undrained friction angle,

φ

u

=0

10 B

z

B=1m

x

y

10 B

(b) Three-dimensional analysis Fig. 2. Problem domain and material properties.

Table 1 Properties of the MEM grids. Grid

Field nodes

Background points

Background cells

Gauss points

2D - Coarse 2D - Fine 3D

1817 14,075 1331

441 1681 1331

400 1600 1000

10,000 78,400 27,000

ME: refers to numerical stress integration using the Modified Euler approach. RM: refers to analytical stress integration using the Return Mapping approach.

5.1. Plane strain analysis In the first instance, a two-dimensional strip footing problem is considered with the properties as represented in Fig. 2(a), with B denoting the width of the strip footing. Due to symmetry, only the right-hand half of the problem has been considered in the analysis. In addition, the footing was subjected to a prescribed overall vertical displacement equal to its total width. According to Prandtl’s plasticity solution [25], the undrained bearing capacity

of a weightless layer of soil under a rigid strip footing, qu , is given by

qu ¼ Nc su ;

ð38Þ

where N c denotes the capacity factor and is equal to 2 þ p  5:14. In order to verify the convergence of MEM with explicit as well as implicit solution schemes, load-displacement curves are plotted in Fig. 3. According to Fig. 3, the MEM method provides robust convergence toward the failure load. In addition, the implicit and explicit solution procedures demonstrate almost identical convergence behaviour. Table 2 shows the predicted capacity factor, the error in calculating the capacity factor, the CPU time normalised by the CPU time of the fastest analysis for the same MEM grid, and the total number of iterations to achieve convergence for each combination of

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

Applied pressure / su

74

6

5.2. Axi-symmetric analysis in 2D

5

The second footing problem considered here is a circular footing through an axi-symmetric analysis. The problem properties are similar to those summarised in Fig. 2(a), with B denoting the diameter of the circular footing, with the total prescribed vertical displacement being equal to B. According to Cox et al. [9], the exact theoretical value for the undrained capacity factor, N c , of soil under a rough circular footing is equal to 6.05. Fig. 4 illustrates the load-displacement relations for the axisymmetric analyses for various solution schemes and MEM grids. According to Fig. 4, the MEM method provides robust convergence toward the failure load. In addition, the implicit and explicit solution procedures again demonstrate almost identical convergence behaviour. Table 3 shows the predicted capacity factor, the error in calculating the capacity factor, the CPU time normalised by the CPU time of the fastest analysis for the same MEM grid, and the total number of iterations to achieve convergence associated with each combination of solution scheme and MEM grid for the axi-symmetric analyses. It can be seen from Table 3 that the bearing capacity of the circular footing is predicted with acceptable accuracy in accordance with the coarseness of the MEM grid. Using the coarse MEM grid, the implicit solver is able to provide a more accurate estimate of the bearing capacity while, in the case of the fine grid, the accuracy of the predicted bearing capacity is comparable for the implicit and explicit solvers. As for the plane strain analyses, the explicit solution scheme is able to predict the bearing capacity factor with almost the same accuracy for different numbers of load increments, while the implicit solver fails to converge for the analyses with a single or 10 load increments. With regard to the computational times, once again it is observed that for equal number of load increments (100), the implicit solution scheme is faster. However, overall, the single-increment explicit solver with analytical stress integration leads to the fastest computational times for both the coarse and the fine grids. Moreover, Table 3 shows that utilising analytical stress integration can reduce the computational times for all solution configurations. This reduction is again more significant for the explicit solver.

4 Grid1 - Implicit Solver (NR) 3

Grid1 - Explicit Solver (ADR) Grid2 - Implicit Solver (NR)

2

Grid2 - Explicit Solver (ADR) Analytical Solution

1 0

0

0.2

0.4

0.6

0.8

1

Displacement / B Fig. 3. Load-displacement curves – plane strain analysis.

solution schemes and MEM grids. The results in Table 2 verify the robustness of MEM with various underlying solution procedures. As expected, the error in the calculation of the bearing capacity decreases by employing a finer grid. It is also seen that for the fine MEM grid, the implicit time integrator can predict the bearing capacity slightly more accurately. An interesting capability of the explicit solver is that it can predict the bearing capacity factor with almost the same accuracy for different numbers of load increments, while the implicit solver fails to converge for analyses with a single or 10 load increments. With regard to the computational time, it is observed that for equal number of load increments (100), the implicit solution scheme is faster. This is due to much fewer iterations being performed in the Newton-Raphson solution compared to the ADR. However, from the point of view of predicting the ultimate bearing capacity, compared to the implicit solver, the single-increment explicit scheme is only slightly more time consuming for the coarse grid while being faster for the fine grid. Moreover, Table 2 shows that utilising the analytical stress integration scheme reduces the computational times for all configurations where a solution was found. This reduction is more significant for the explicit solver due to the larger number of iterations that it needs.

Table 2 MEM results for plane strain analyses. Grid

Solution scheme

Load increments

Nc

Error in N c (%)

Total CPU time (normalised)

Total iterations

Coarse

NR-ME

1 10 100 1 10 100 1 10 100 1 10 100

NCa NCa 5.31 NCa NCa 5.31 5.31 5.31 5.31 5.31 5.31 5.31

– – 3.31 – – 3.31 3.31 3.31 3.31 3.31 3.31 3.31

– – 1.00 – – 1.00 1.84 6.25 14.78 1.64 2.47 7.41

– – 262 – – 270 9250 11,992 42,344 9066 12,067 42,537

1 10 100 1 10 100 1 10 100 1 10 100

NCa NCa 5.14 NCa NCa 5.14 5.14 5.15 5.15 5.14 5.15 5.15

– – 0.00 – – 0.00 0.00 0.19 0.19 0.00 0.19 0.19

– – 1.27 – – 1.24 1.21 4.02 7.03 1.00 1.94 4.21

– – 308 – – 315 20,317 27,012 69,416 23,450 33,818 78,536

NR-RM

ADR-ME

ADR-RM

Fine

NR-ME

NR-RM

ADR-ME

ADR-RM

a

NC: iterations not converged to the failure load using the specified number of load increments.

75

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

6

5

Applied pressure / su

Applied pressure / su

6

Grid1 - Implicit Solver (NR)

4

Grid1 - Explicit Solver (ADR)

3

Grid2 - Implicit Solver (NR) Grid2 - Explicit Solver (ADR)

2

Analytical Solution 1

5 4 Implicit Solver (NR)

3

Explicit Solver (ADR) 2

Analytical Solution

1

0

0

0.2

0.4

0.6

0.8

1

0

Displacement / B

0

0.2

0.4

0.6

0.8

1

Displacement / B

Fig. 4. Load-displacement curves – axi-symmetric analysis. Fig. 5. Load-displacement curves – three-dimensional analysis.

5.3. 3D strip footing The final example involves three-dimensional analysis of the undrained bearing capacity of soil under a rigid rough strip footing subjected to vertical prescribed displacement. The properties of the problem are summarised in Fig. 2(b) and the vertical prescribed displacement is set equal to the total width of the footing. Fig. 5 shows the load-displacement curves for the three-dimensional analyses with the implicit and explicit solution procedures. According to this figure, the MEM method provides robust convergence toward the failure load. As in the previous examples, the implicit and explicit solution procedures demonstrate almost identical convergence behaviour. Table 4 presents the predicted capacity factor, the error in calculating the capacity factor based on Prandtl’s exact solution, the CPU time normalised by the CPU time of the fastest analysis for the same MEM grid, and the total number of iterations to achieve convergence associated with each combination of solution schemes and the MEM. This Table reveals that the bearing capacity is predicted by the explicit and implicit solvers with almost the same accuracy. The explicit solution scheme is again more robust and is able to predict the bearing capacity factor with almost the same accuracy for different numbers of load increments, while

the implicit solver fails to converge the analyses with less than or equal to 10 load increments. The single-increment explicit solver with analytical stress integration gives the fastest computational time, and the analytical stress integration scheme reduces the computational times for all configurations where a solution was obtained. This reduction is more significant for the explicit solver due to its larger number of iterations. 5.4. Large scale computational efficiency As seen in the numerical results from three footing sample problems, while the explicit time integration scheme is capable of estimating the failure load in only a single load increment, for 100 load increments it is much slower than the implicit time integration scheme. However, one important observation about the speed of the implicit and explicit solution approaches, which is not immediately obvious from the numerical results shown, is that the computation times per iteration grows significantly for the implicit method as the problem sizes increases, while such growth is minimal for the explicit approach. Fig. 6 illustrates the relation between the size of the problem and the CPU time per iteration for various solution scheme configurations. Note that the total number of shape function values for all the integration points is

Table 3 MEM results for axi-symmetric analyses. Grid

Solution scheme

Load increments

Nc

Error in N c (%)

Normalised CPU time

Total iterations

Coarse

NR-ME

1 10 100 1 10 100 1 10 100 1 10 100

NCa NCa 6.22 NCa NCa 6.22 6.26 6.27 6.27 6.26 6.27 6.27

– – 2.81 – – 2.81 3.47 3.64 3.64 3.47 3.64 3.64

– – 1.12 – – 1.11 1.08 7.26 14.54 1.00 2.36 10.42

– – 241 – – 249 3105 7069 31,076 2968 7107 31,207

1 10 100 1 10 100 1 10 100 1 10 100

NCa NCa 6.08 NCa NCa 6.08 6.09 6.1 6.1 6.09 6.1 6.1

– – 0.50 – – 0.50 0.66 0.83 0.83 0.66 0.83 0.83

– – 3.12 – – 3.09 1.48 5.70 15.29 1.00 2.41 7.61

– – 262 – – 270 5843 14,516 49,730 6464 15,675 49,868

NR-RM

ADR-ME

ADR-RM

Fine

NR-ME

NR-RM

ADR-ME

ADR-RM

a

NC: iterations not converged to the failure load using the specified number of load increments.

76

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77

Table 4 MEM results for 3D analysis. Solution scheme

a

Load increments

Nc a

Error in N c (%)

Normalised CPU time

Total iterations

NR-ME

1 10 100

NC NCa 5.53

– – 7.59

– – 6.43

– – 243

NR-RM

1 10 100

NCa NCa 5.53

– – 7.59

– – 6.04

– – 267

ADR-ME

1 10 100

5.51 5.52 5.53

7.20 7.39 7.59

1.31 9.80 23.11

5747 15,928 57,220

ADR-RM

1 10 100

5.51 5.53 5.54

7.20 7.59 7.78

1.00 4.04 13.45

6329 16,432 60,875

NC: iterations not converged to the failure load using the specified number of load increments.

CPU time per iteration (s)

80

NR-ME NR-RM ADR-ME

60

ADR-RM

40

20

0 0

10

20

30

40

50

Normalised problem size Fig. 6. CPU time per iteration for growing problem sizes.

utilised as an indicator of the problem size in this figure. Other potential indicators of the size of the problem include the bandwidth of the global stiffness matrix or the maximum size of the local stiffness matrices over all integration points. According to Fig. 6, the CPU time for each iteration of the implicit solver grows much more significantly than the CPU time of the explicit approach. This implies that, for larger problems, the implicit solver will tend to exhibit comparable computational times to those of the explicit solver. In particular, for three-dimensional problems, the cost of factorizing and solving the resulting linear system of equations in each iteration of the implicit solver increases dramatically. 6. Conclusion In this paper, the application of the Maximum Entropy Meshfree (MEM) method to elastoplastic analysis of geotechnical problems was studied. Explicit and implicit approaches were investigated for solution of governing equations with either numerical or analytical schemes for stress integration. The robustness and stability of the MEM method, along with four configurations of the solution procedures, were verified by performing two- and three-dimensional numerical experiments, including the bearing capacity of soil under rigid and rough footings. Overall, the following remarks could be concluded from this study:  For the nonlinear geomechanics problems considered in this study, the MEM method provides a stable and robust numerical tool. All solution configurations provide stable convergence results.

 The explicit time integration scheme allows prediction of the failure load for analyses with any number of load increments, including only a single load increment, while the implicit approach fails to converge to the failure load with large load increments.  Based on our numerical results, where both the explicit and implicit approaches converge to the failure load, the implicit solver is faster in terms of CPU time.  It is observed that the CPU time per iteration grows rapidly for the implicit approach, while such an increase is minimal for the explicit approach. This implies that, for larger problems, the explicit solver may outperform the implicit solver in terms of total computational time for the same number of load increments.  The analytical stress integration scheme was shown to be robust and to accelerate the MEM analysis. The acceleration tends to be more substantial in case of the explicit scheme since a large number of iterations is involved.

References [1] Abbo AJ. Finite element algorithms for elastoplasticity and consolidation PhD Thesis. Australia: The University of Newcastle; 1997. [2] Alamatian J. A new formulation for fictitious mass of the dynamic relaxation method with kinetic damping. Comput Struct 2012;90:42–54. [3] Arroyo M, Ortiz M. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods. Int J Numer Meth Eng 2006;65(13):2167–202. [4] Bathe KJ. Finite element procedures. Klaus-Jurgen Bathe; 2006. [5] Beissel S, Belytschko T. Nodal integration of the Element-Free Galerkin method. Comput Methods Appl Mech Eng 1996;139:49–74. [6] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Meth Eng 1994;3:229–56. [7] Belytschko T, Mullen R. Explicit integration of structural problems. In: Bergan P et al., editors. Finite elements in nonlinear solid and structural mechanics, Trondheim, Norway; 1977. [8] Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics. IBM J Res Dev 1967;11(2):215–34. [9] Cox AD, Eason G, Hopkins HG. Axially symmetric plastic deformations in soils. Philosoph Transact Roy Soc Lond A: Math, Phys Eng Sci 1961;254(1036):1–45. [10] Dang HK, Meguid MA. Evaluating the performance of an explicit dynamic relaxation technique in analyzing non-linear geotechnical engineering problems. Comput Geotech 2010;37(1):125–31. [11] Jaynes ET. Information theory and statistical mechanics. Phys Rev 1957;106:620–30. [12] Kardani O, Krabbenhøft K, Lyamin AV. Application of adaptive dynamic relaxation to highly nonlinear geotechnical problems. In: 11th World Congress on Computational Mechanics (WCCM XI), Barcelona, Spain, July 2014. [13] Khoshghalb A, Khalili N. A meshfree method for fully coupled analysis of flow and deformation in unsaturated porous media. Int J Numer Anal Meth Geomech 2013;37(7):716–43. [14] Kreig RD, Key S. Transient shell response by numerical time integration. Int J Numer Methods Eng 1973;17:273–86. [15] Krieg RD. Unconditional stability in numerical time integration methods. J Appl Mech 1973;40(2):417–21. [16] Liu GR, Gu YT. An introduction to meshfree methods and their programming. The Netherlands: Springer; 2005.

O. Kardani et al. / Computers and Geotechnics 84 (2017) 68–77 [17] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. Int J Numer Meth Eng 2001;50(4):937–51. [18] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246(1):29–46. [19] Millan D, Rosolen A, Arroyo M. Thin shell analysis from scattered points with maximum-entropy approximants. Int J Numer Meth Eng 2011;85(6):723–51. [20] Nazem M, Kardani M, Bienen B, Cassidy M. A stable Maximum-Entropy Meshless method for analysis of porous media. Comput Geotech 2016;80:248–60. [21] Oakley DR, Knight NF. Adaptive dynamic relaxation algorithm for non-linear hyperelastic structures Part I. Formulation. Comput Methods Appl Mech Eng 1995;126(1):67–89. [22] Oliaei MN, Soga K, Pak A. Some numerical issues using element-free Galerkin mesh-less method for coupled hydro-mechanical problems. Int J Numer Anal Meth Geomech 2009;33:915–38. [23] Ortiz A, Puso M, Sukumar N. Maximum-entropy meshfree method for incompressible media problems. Finite Elem Anal Des 2011;47(6):572–85. [24] Park KC, Underwood PG. A variable-step central difference method for structural dynamics analysis—Part 1. Theoretical aspects. Comput Methods Appl Mech Eng 1980;22(2):241–58. [25] Prandtl L. Hauptaufsätze: Über die eindringungsfestigkeit (härte) plastischer baustoffe und die festigkeit von schneiden. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 1921;1(1):15–20. [26] Quaranta G, Kunnath S, Sukumar N. Maximum-entropy meshfree method for nonlinear static analysis of planar reinforced concrete structures. Eng Struct 2012;42:179–89.

77

[27] Rezaiee Pajand M, Alamatian J. The dynamic relaxation method using new formulation for fictitious mass and damping. Struct Eng Mech 2010;34:109–33. [28] Shibata T, Murakami A. A stabilization procedure for soil-water coupled problems using the element-free Galerkin method. Comput Geotech 2011;38:585–97. [29] Sloan SW, Abbo AJ, Sheng DC. Refined explicit integration of elastoplastic models with automatic error control. Eng Comput 2001;18(1/2):121–54. Erratum: Engineering Computations, 19(5/6), 594–594, 2002. [30] Sukumar N. Construction of polygonal interpolants: a maximum entropy approach. Int J Numer Meth Eng 2004;61:2159–81. [31] Sukumar N, Wright WR. Overview and construction of meshfree basis functions: from moving least squares to entropy approximation. Int J Numer Meth Eng 2007;70:181–205. [32] Ullah Z, Coombs WM, Augarde CE. An adaptive finite element/meshless coupled method based on local maximum entropy shape functions for linear and nonlinear problems. Comput Methods Appl Mech Eng 2013;267:11–132. [33] Underwood P. Dynamic relaxation. In: Belytschko T, Hughes TJR, editors. Computational methods for transient dynamic analysis. North Holland, Amsterdam; 1983. p. 246–65. [34] Wang JG, Karim MR, Lin PZ. Analysis of seabed instability using element free Galerkin method. Ocean Eng 2007;34:247–60. [35] Wang JG, Liu GR, Lin P. Numerical analysis of Biot’s consolidation process by radial point interpolation method. Int J Solids Struct 2002;39(6):1557–73. [36] Zienkiewicz OC, Valliappan S, King IP. Elasto-plastic solutions of engineering problems ‘initial stress’, finite element approach. Int J Numer Meth Eng 1969;1 (1):75–100.