Arbitrary Lagrangian-Eulerian formulation for element-free Galerkin method

__k!!B

Ck-QA ELSEYIER

Computer methods in applied mechanics and englneerlng Comput.

Methods

Appl. Mech. Engrg.

152 (1998)

19-46

Arbitrary Lagrangian-Eulerian formulation for element-free Galerkin method J.-P. Ponthot*, Department

T. Belytschko

of Civil and Mechanical Engineering, Robert McCormick School of Engineering and Applied Sciences, Institute, Northwestern University, Evanston, IL 60208-3109, USA Received

26 February

The Technological

1997

Abstract Arbitrary Lagrangian-Eulerian (ALE) formulation of the Element Free Galerkin (EFG) method is presented. EFG is a meshless method for solving partial differential equations in which the trial and test functions employed in the discretization process result from moving least square interpolants. The most significant advantage of the method is that it requires only nodes and a description of internal and external boundaries and interfaces, such as cracks, of the model: no element connectivity is needed. However, as for any discretization method, acceptable solutions can only be obtained for a sufficiently refined discretization. In dynamic fracture problems, where the crack path can be arbitrary, and is thus a priori unknown, this necessitates a refined discretization in large parts of the computational domain which can lead to prohibitive computation ccIsts. ALE formulation allows to continuously relocate nodes on the computational domain. By combining EFG with ALE, it is thus possibie, in a crack propagation problem, to refine locally the spatial discretization in the neighborhood of a propagating crack-tip. Results are presented for a wave propagation problem as well as for 2-D dynamic crack propagation problems.

1. Introduction For the past 25 years or so, the finite element method (FEM) has been an industrial standard tool for solving a wide variety of mechanics problems. However, for some highly nonlinear problems where the geometry can evolve significantly such as large deformation problems or fracture dynamics problems with propagating cracks, some difficulties are not yet overcome. These difficulties generally arise due to the inherent structure of the FEM and topological knowledge it requires, namely the rigid connectivity defined by elements. In fracture problems, for instance, finite element edges provide natural lines along which cracks can grow. However, if the crack path is not known a priori, FEM requires remeshing in order to follow an arbitrary crack path. In contrast to the FEM, meshless methods do not have a rigid connectivity provided a priori, making them ideal in applications where finite element methods have the most difficulties. One meshless method, the element free Galerkin (EFG) method, proposed by Belytschko and coworkers [8,12-171 and based on moving least square (MLS) approximations, has been used extensively to solve fracture problems where it offers considerable simplifications. In EFG, the connectivity and shape functions are constructed by the method and they can be automatically modified during the computation to account for crack propagation. There are several meshless methods: Smoothed Particle Hydrodynamics-Gingold and Monaghan [33,68], Generalized Finite Ljifference Method-Liszka and Orkisz [57], Diffuse Element Method-Nayroles et al. [72-751, Wavelet Galerkin Method-Qian and Weiss [86], Multiquadrics-Kansa [47,48], Reproducing Kernel Particle Method-Lin et al. [58-601, Element Free Galerkin-Belytschko et al. [ 12,14,17] and Hegen [36], the

* Corresponding

author.

0045-7825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PZI SOO45-7825(97)00180-l

20

J.-P. Ponthot, T. Belytschko

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

Partition of Unity FEM or PUFEM--BabuSka and Melenk [5,6], Hp Clouds-Duarte and Oden [28,27] and Finite Point Method-Ofiate et al. [79]. However, most of them make use of a weight function of compact support and many of them, explicitly or implicitly use MLS approximations. From the Continuum Mechanics viewpoint, two classical formulations co-exist in the literature: in solid mechanics, the Lagrangian description is employed extensively to describe the evolution of a given set of material particles whereas the Eulerian description is preferred when it is convenient to model a fixed region in space, as is generally the case in fluid mechanics. However, it has long been recognized that for many engineering problems, such as fluid-structure interactions, impact, contact, penetration and metal-forming simulation, neither the Lagrangian nor the Eulerian descriptions are well suited because of their well-known intrinsic limitations. Essentially, the Eulerian description can handle very large material distortions but have limited spatial resolution to track a moving boundary while the Lagrangian description can precisely track a moving boundary or interface but is handicapped by mesh distortions. In fact, the Eulerian and Lagrangian formulations viewpoints have rather complementary virtues. Therefore, some researchers have tried to establish a new formulation called Arbitrary Lagrangian-Eulerian (ALE) that combines best of both. In 1964 were published the first two papers dealing with ALE. Noh [76] and Frank and Lazarus [30] made contributions in that sense, for fluid mechanics problems, under the name of Coupled Eulerian-Lagrangian and Mixed Eulerian-Lagrangian, respectively. They implemented their formulations in a finite difference code. The finite element community originally adopted ALE formulation for fluid-structure interaction calculations problems [3,7,10,11,24-26,34,43,49,61,65,77,78]. It was further extended to solid mechanics by various researchers among which: [1,2,4,18,19,22,32,35,37-40,51,52,62-64,81-85,87,91]. In this paper, we present an EFG formulation in an ALE framework. This formulation combines both the advantages of EFG and ALE over more classical FEM/Lagrangian description. One of the most attractive features of EFG is that moving cracks can be simply modeled by an extending line. Examples of such simulations in a Lagrangian framework can be found in [ 12-14,171. By contrast, FEM encounters considerable difficulties in treating a moving crack. However, as is the case for any discretization method, a more precise solution can only be achieved in EFG by increasing the nodal density where appropriate. In previous EFG simulations of dynamically growing cracks, where the path is a priori unknown, a finer discretization was used on the whole region of interest, thus leading to much higher computational costs. The ALE formulation is an attractive alternative in the sense that, as will be shown hereafter, it is possible to superpose on a rather coarse Lagrangian discretization a limited number of ALE nodes located around the crack that move at the same velocity as the crack tip. Thus, the density of nodes can be kept dynamically high in a limited, though evolving, region where it is needed and at a rather small cost (few nodes are involved in the ALE process). Moreover, the motion of the ALE nodes, follows directly from the crack propagation and thus can be easily and automatically determined without recourse to any more or less elaborate mesh management procedure. The user intervention is limited to the pre-processing phase where he has to generate the initial position of the ALE nodes.

1. Element-free

Galerkin

approximation

For completeness, we shall recall here the fundamentals of the Element-Free Galerkin method. More details about this method can be found in the work of Belytschko and coworkers [8,12-171. The EFG method is a meshless or particle method. The objective of such a meshless method is to obtain an approximation to a function u(x, t) strictly in terms of parameters at a set of nodes and a description of the boundaries of the domain of interest. This boundary description includes any interior surfaces of discontinuity such as cracks or interfaces between different materials. In the EFG method, the field variable ~(1, t) is approximated by MLS approximations [55]. The MLS approximations u~(x, t) is given by

uh(x,t) =

i:

p,(x)a,(x, t)

=pT(x>a(x, t)

,=I

where p(x) is a complete polynomial

basis of arbitrary order and a(x, t) are coefficients

(to be determined)

which

J.-P. Ponthot, T. Belytschko

are functions of the space coordinate following basis: PT(N = 11

x

y]

pT@) = [l

x

y

x and time t. For example,

in two dimensions,

one can choose

(m = 3, linear basis) xy

x2

y*]

21

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

the

(2)

(m = 6, quadratic basis)

(3)

Actually, p(x) is not restricted to be a polynomial and the basis can be enhanced by non-polynomial terms. Examples using such an enriched basis can be found in [29]. The unknown coefficients a&, t) in (1) are obtained at any point x by minimizing the following weighted, discrete error norm (4) /=I

where w(x - x,) is a weight function of compact support (often called the domain of influence of node Z) and R is the number of nodes whose support includes point x. The elements of the discrete set x,, I = 1, . . . , n of nodes such that w(x -x,) 3’0 are called the neighbor nodes at point x (see Fig. 1). In expression (4), u,(t) is the parameter associated with node I of the approximated field; it is not the nodal value because u”(x, t) is an approximation, not an interpolation. Minimization of (4) with respect to a(x, t) then yields to the following system of linear equations for the coefficient Q(X, t): A(x)u(x,

t) = B(x)u(t)

(5)

A(x) = ,g, w(x - x,)P(x,)P~(x~)

(6)

where

B(x) = [w(x -x,)p(x,) UT(l) = [u,(t)

$0)

w(x -x&(x*)

...

w(x -x,,)p(x,,)I

(7)

(8)

. . . u,,(t)1

where u(t) is the vector of nodal unknowns. If A is invertible (a necessary but not sufficient condition is that n 2 m; i.e. the number of neighbors is at least equal to the basis dimension), the coefficients a(x, t) can be expressed as u(x, t) =A-‘(x:B(x)u(t) By substituting

(9)

(9) into (l), the MLS approximates

are finally given by

NODE /

EVALUATION /

POINT

Fig. I. Illustration of nodal domains of influence in two dimensions. The neighbor list for evaluation since their domain of influence contains point x. Node 3 is excluded from the neighbor list for x.

at point x includes nodes

1, 2 and 4

22

J.-P. Ponthot, T. Belytschko

L&X, t) =pT(x)a(x,

/ Compuz. Methods Appl. Me&

Engrg. 152 (1998) 19-46

t) = GT(X)U(l)

(10)

or

The EFG shape function

@(x, t) is defined as

(12)

GT(x) =pT(x)[A-‘(x)B(x)]

Thus, in EFG, the list of neighbors and the resulting shape functions can be automatically evaluated once a given weight function w has been chosen. The continuity of the shape functions is governed by the continuity of the basis function p(x) and the smoothness of the matrices A-‘(x) and B(x). The latter is governed by the smoothness of the weight functions (see [54,66,74]). It should also be noted that even if W(X- x,) is a polynomial, the shape functions are no longer polynomials and that the shape functions do not satisfy the Kronecker delta condition at each node, i.e. (13)

Q/(x,) f a,,

Consequently, the imposition of essential boundary conditions is more complicated than for the standard FEM. Several methods have been developed by Belytschko and coworkers, including Lagrange multipliers [8,15], penalty method [ 131, modified variational principles [661, collocation at nodes for explicit dynamics [ 171, weak form of the kinematic boundary conditions [67,14] and FEM-EFG coupling [ 16,53,80]. Moreover, in a crack propagation problem, the weight functions w(x -x,) at a given point are functions of time because of the support of a node change with time in the vicinity of the crack. To account for that effect, the shape functions have to be periodically reevaluated, at least in the vicinity of a crack tip. I. 1. Weight functions In the absence

of a crack, the weight function

associated

with the Ith node is non-zero

B, = {x E R* : Ix -x,1 s r,} where r, is the norm. The ball support of @,. We define a r_

(14)

radius of support of the Zth node, and Ix -x,1 is the distance between x and x, in any suitable B, defines the largest set on which the shape functions @, may be non-zero, i.e. it defines the In other words, IA, influences the approximation of u only on B,. normalized distance:

lb-4

(15)

YI and require the weight function w(r)

30 =0

1

w(r) to be non-negative

inside the ball B,, i.e.

for r < 1 forral

(16)

For EFG, the weight function w(r) is generally chosen as a monotonically distance r. Here, we have adopted the following Gaussian distribution: -w/0.4?

w(r) = e 10

in a ball

for

r

<

decreasing

function of the normalized

1

for r 2 1

(17)

Other examples of weight functions including cubic splines and polynomials can be found in [16,80]. Those weight functions are all defined as functions of the single parameter r. If the ray from x to xr does not intersect any boundaries, the simplest definition of r is given by the normalized Euclidean distance between two points, i.e.

23

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

(a)

concaveCUrYe

(b)Crackor

corner

Fig. 2. Scheme for computing the weight function parameter r.

(18) However, for non-convex domains, it is necessary to modify the standard definition of r to account for non-convex boundaries (cracks, concave curve, corner, . . .). The simplest approach is to cut off the weight functions by the boundary in an abrupt manner by setting the weight to zero at points not visible from the node (visibiliry crirerion). The so defined weight function and thus the resulting shape functions are no longer Co. In order to overcome that drawback, Belytschko and coworkers [12,80] have redefined parameter r by using a measure of distance along a path which lies entirely within the domain, illustrated in Fig. 2, which we have used. This is called the diffraction method. In the diffraction method, the expression for Y is replaced by an expression for the distance along a path consisting of two straight lines which lie entirely within the domain so that:

y= 1%- TII + Ik -xcII rl

(19)

where xc is a point selected so that the two rectilinear line segments lie entirely within the domain. In this case, the support of the shape functions is defined, not by (14), but by: B, = {x E R* : w(x, x,) > 0) which allows the weight function support to be of almost any shape. Krysl and Belytschko [54] have discussed the convergence and conformity (diffraction method) and those by the visibility criterion.

(20)

for these

shape

functions

2. ALE basic equations 2.1. Notations

and preliminaries

The Lagrangian or material coordinate X is considered as a ‘marker’, attached to a given material particle. It describes its motion through a function q which gives the position of the particle as a function of time t: x = 44-F t)

(21)

where x is the spatial coordinate (the current position of the material particle identified by X) and +J is the material motion function or the mapping function between X and x at a given time t. Alternatively, it is also possible to describe the motion of a continuum by a set of coordinates X, called the

24

J.-P. Ponthot, T. Belytschko I Comput.Methods Appl. Mech. Engrg. 152 (1998) 19-46

referential coordinates. Generally, the motion of the referential coordinate point can be arbitrary and therefore not coincide with the material coordinates. With this set of coordinates, the motion of the continuum can be described by a relation analogous to (21) i.e. x = ax,

t)

(22)

where @ is the referential motion or the mapping between x and X. The material and referential and ti are respectively defined as u=x-x;

displacements

ti=x-X

u

(23)

Though the referential coordinates are arbitrary, they have to be regular in the sense that, at every time t, a one-to-one (though not always the same) relation must exist between the current configuration 0 of the body and the set of referential coordinates fix. This is guaranteed by .?=det@#O

VXEQ~

where

This relation is totally analogous configuration Q,: J=detF#O

VXXE,,

g=$

(24)

to the classical relation that holds in a Lagrangian

where

on the reference

F=$

(25)

Relations (24) and (25) imply that, at each time t, there is also a one-to-one coordinates and the material coordinates. 2.2. Time derivative

description

relation

between

the referential

of a function

Generally, to describe the behavior of a continuum, time derivative of some function f is required. As any function f can be described as a function of x, X or x and three types of time derivatives can be defined, namely: The spatial derivative df V&9 t) -_=at at

(26)

x

The material derivative f=

(;lf(x, t) at

The referential

derivative

p= Jf(X> t) at For example,

(27)

x

(28)

x

the material velocity is classically

defined by (29)

whereas, the referential

velocity can be defined by (30)

Those time derivatives are not independent. The fundamental relationship derivative can be established readily, for any quantity f as [43,25,84]:

f=fp+c.vf=p+ci-g I

between material and referential

time

(31)

25

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

where c is the relarive

or convective

velocity between

the material

u and referential

ri velocity:

c=v-6

(32)

Actually, the term c *Vf in (3 1) represents a convective effect due to the relative motion between the material and referential descriptions. As D is arbitrary, the ALE formulation has generally no intrinsic dependence on particles or position. However, the following basic viewpoints may be stressed: (1) B = 0: the reference frame is fixed in space and this corresponds to the Eulerian viewpoint. (2) 6 = u: the reference frame moves (locally) in space at the same velocity as the particles and this corresponds to the Lagrangian viewpoint. (3) 0 # ti # u: the reference frame moves in space at an arbitrary velocity. This is the most general case termed Arbitrary Lagrangian-Eulerian. 2.3. Governing

equations

of the continuum

We consider the transient problem in elasticity with small displacements the material description, the equations of motion are q,,, + pb, = pti,

on the domain

in R

0 bounded by X In

(33)

where (T is the Cauchy stress tensor corresponding to the displacement field U; ti is the material acceleration; b is the body force velstor and p is the mass density. We use the usual tensor notation in which repeated indices are summed over the space dimension and, in (33), an index preceded by a comma denotes the derivative with respect to spatial coordinates. The corresponding referential equation is readily obtained from (31) and (33). It reads aij,j + pb, = ~(6~ + c~u~,~) In (33) and (34), the Cauchy

in 0

stress tensor is obtained

(34) from the elastic constitutive

law

q,i = H;jklEkl

(35)

and the strain-displacement ‘kl

=

;

@kJ

+

equation

U,,k)

The weak form corresponding to the conservation equation is obtained by multiplying admissible test function SU over the spatial domain 0 and employing the divergence traction forces on the boundary 2. This results in

(36) (34) by kinematically theorem to imbibe the

where t, = gjjni are the natural boundary conditions applied on 2,. The discretized system of equations is obtained by substituting the following expressions for, respectively, material displacement, test functions, convective velocity, material velocity, referential acceleration and initial velocity, into (37): (38)

(39)

c(xt t) = v(x

2 @/(X>C,@>

t) = c

I=1

45(x)u,(t)

(40)

(41)

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

26

q(x>a,@)

d(x, t) = i

(42)

,=I

u,(x) = i

,=I

qs(x)u”,

(43)

where @, are the EFG shape functions given by ( 12) where the spatial coordinates x have been replaced by the referential coordinates x. The assumption made in (24) guarantees that this change of variable is regular. REMARK 1. The shape functions are not computed in terms of referential operator-split methodology is used (see Section 3.1). Introduction

of (38)-(43)

M; + CV +fint

coordinates

in the program because an

into (37) results in

zf”“’

(44)

with f In’= Ku and W, = i I, p@,@,I dR

(45)

C,, = I PC, a: @,:.,I da n

(46)

K,, = I

B TDB, da

(47)

f;‘=i,B;udD

(48)

f 7’ = j-, p@,b dR + i @,t d& 9,

(49)

where M is the mass matrix, C the convective matrix, K the stiffness matrix, f I”’ are the internal forces and f “’ are the applied forces. In two dimensions, we have

(50)

(51) and, for plane strain

pu E D=(1+V)(l-2V)

u

v o

l-v OL

where E and Y are the Young’s modulus

0 0 I-2u

1

2

-l

(52)

I

and Poisson ratio, respectively,

and Z is the 2 X 2 identity

matrix.

REMARK 2. In dynamic fracture problems, as the crack propagates, the weight functions are functions of time because the domain of influence changes. As a consequence, the shape functions are also time dependent. Thus, the separation of variables, which was implicitly supposed in (38)-(43) is not strictly valid. However, the contribution of the shape functions time derivatives has been neglected in the present approach. Nevertheless, at each time step, when it is necessary, the shape functions are updated and a new mass-matrix is recomputed. This remark holds for both Lagrangian and ALE formulations.

J.-P. Ponthot, T. Belytschko

3. Numerical

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

21

integ,ration

3. I. Operator-split

methodology

One of the major difficulties associated with the time integration of (44) results from the presence of the non-symmetric convective term Cu (meaning the problem is not self-adjoint). This term stems from the relative motion between the material and referential system. In a Lagrangian formulation, this relative motion vanishes, since the material and referential systems are identical. Therefore, no convection occurs and the C matrix is identically zero. In fact, two ways are possible in implementing and solving the ALE equations. These two ways also correspond to the ‘two approaches taken in implementing an Eulerian solver from the viewpoint of fluid mechanics. The first approach solves the fully coupled equations as they are written. In the Finite Element community, Belytschko et al. [7,49,10,11], Hughes et al. [43] and Liu [61,65] have adopted this approach for fluid-structure interactions. For solid mechanics modeling Liu et al. [62-641 have also successfully adopted and extended Streamline Upwind Petrov-Galerkin (SUPG) methodology, as developed by Hughes and coworkers [21,44]. Other examples are given in [87] and [32]. Though all these authors use a coupled formulation for determining u and fi, it must be pointed out that, generally, when an iterative, implicit algorithm is used to solve the equations, the mesh management procedure (i.e. finding a value for the 6’s) is either rather simple (degrees of freedom are a priori declared either Lagrangian or Eulerian [43]), or it is treated outside the iterative loop [64,32]. The only noticeable exception is given by Schreurs et al. [87] under quasi-static assumption. They solve the fully coupled system of equations (i.e. where both u and ti are unknowns) thanks to a fictitious elastic material which is used to govern mesh displacements in order to minimize mesh distortions. Alternatively, one can use an operator splitting to integrate the equations in time. We have adopted this technique here. In fact, this technique is a very convenient method for breaking complex problems into a series of less complicated ones (see e.g. [23] for more details). Typically, a Lagrangian step is performed first. During this phase, the grid moves with the material during the step and its configuration at the end of the step is different from a (to be defined) reference grid. The solution obtained from the Lagrangian phase is then mapped (convected) to the reference grid in order to complete the ALE step. This algorithm can be readily extended to EFG where, instead of moving the mesh, nodes (all of them or some of them) can be moved arbitrarily during the convection step. The main advantage of this operator split is that a complex problem is broken into much simpler ones and that for each of these time substeps, efficient algorithms can be applied. As a consequence, algorithms that rely on operator split are usually robust and less expensive than the algorithms devised for solving the fully coupled equations. In the finite element community, operator split algorithms have been used by Donea [24-261 flor fluid-structure interaction problems as well as by Benson [ 18,191, H&ink and coworkers [40,39,2,1], Huerta et al. [37,38], Baaijens [4] and Ponthot [81-851. Fracture mechanics problems were also addressed by Koh et al. [35,52,51]. Other advantages of the split algorithm will be discussed further.

3.2. Lagrangian

su,bstep

The first substep is rather classical since the material and referential coordinate systems are supposed to be identical. As a consequence, the convective velocity c is trivially zero and material and referential time derivative are identical. Another consequence is that the shape functions Q,(x) are identical to classical Lagrangian shape functions @,(X). Therefore, the discretized equation of momentum conservation classically reads Ma + Ku

=f’“’

where a = ti and a typical time-step now consists of finding the displacement u which satisfies the semi-discrete equation (53) associated with initial conditions U” and un that prevail at the beginning of the current time-step, i.e. at time tn. We have used the explicit central difference algorithm (see e.g. [9] for details) for which we have (the central difference procedure classically includes some asynchronization on the initial condition on the velocity):

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

28 a”

u u

(&,p-‘(f’“’

= ,,+I/*

=U I!+

I

=

uli

-Ku”)

“-‘I2 + Ata” +

(55)

Atu,,+1/2

where At is the time increment such that t”+’ = t” + At and ML is a lumped mass matrix obtained from the consistent one by a row-sum technique. In order to perform the spatial integration for M, K, f’“’ and f ex' m (53), the quadrature scheme with a background cell structure proposed in [ 15,661 is used. The cells are independent of the nodes and are arranged in a regular pattern in all dimensions. In forming the discrete equations, each cell is considered in turn and Gaussian quadrature is performed. Each quadrature point contributes non-zero entries to the equations only to those nodes whose domain of influence contains xo, the quadrature point position or, in other words, only to these nodes I such that w(xo -x,) > 0 or @,(xo) > 0. Thus, in the final equations, no coupling occurs between nodes which are not within a specified distance (which corresponds to the domain of influence) and the corresponding entries in K and M vanish. Since the mass matrix M is symmetric, positive-definite and K is symmetric, positive semi-definite, the consistency, convergence and stability analysis of time integration procedures for finite element methods (see [9] for an account) also hold in the implementation of the EFG/ALE method. Therefore, the critical time-step that can be used in this part of the algorithm is totally independent of the amount of convection that takes place in the second part of the algorithm. This is certainly not the case of the fully coupled algorithm where the critical time-step that can be used also depends on the convective phase, as will be illustrated in the wave propagation problem. Another advantage of the operator split is that any available algorithm e.g. Newmark, Hilber-Hughes-Taylor or any other [41], can be substituted to the explicit central difference scheme described here without affecting much the Eulerian part of the algorithm. Thus, the operator split methodology is very versatile. A third advantage is that robust and well-suited algorithms developed for Lagrangian codes can be used without any modification, thus retaining their good accuracy and effectiveness. For example in the evaluation of the stress intensity factor for fracture mechanics problems, the J-integral can be evaluated by using reliable domain integrals as described in [69,70]. Another example is given by the large deformation/displacement domain where objective constitutive equations can be time integrated by powerful schemes that retain the objectivity property for arbitrary large time-steps (incremental objectivity) (see e.g. [42,84]). If the incremental objectivity property is generally not thought to be an issue for explicit codes where the time-steps are very small (however, a good counter-example in the large displacement/small strain domain can be found in [50]), it is fundamental for implicit algorithms where the time-steps are rather large. In such a case, if the incremental objectivity property is not achieved, important spurious stresses can result. The disadvantage of the operator split methodology is that, as shown in [23,19,88], the splitting results in a first-order accurate solution, even if both Lagrangian and Eulerian phases are second-order accurate. However, as pointed out by Benson [18], this argument may only be theoretical since the Lagrangian step is rarely fully second-order accurate in time. For example, the radial return procedure [g&42,84] which is often used to integrate elasto-plastic constitutive equations is only first order accurate. It should also be noted that the fully coupled explicit algorithm also includes some inherent asynchronization [7] since the transport term evaluation in (44) is based on V”-I’* rather than u”. Therefore, second-order accuracy cannot be guaranteed in this case either. Moreover, arbitrarily moving some nodes with a non-zero convective velocity reintroduces some time-dependency of the weight and shape functions (see Remark 2). As expressions for the shape function time derivatives are not available, this should be considered as another source of error. 3.3. Eulerian

substep

Once a new equilibrated configuration has been reached from the Lagrangian takes place. First of all, a new nodal pattern has to be defined. This new nodal (the simplest case is the Eulerian formulation, the displacement associated to the be automatically recomputed in order to achieve some goal. In EFG, we will try to achieve and retain some higher node concentration in

substep, the Eulerian substep pattern can be user-predefined nodes is always zero) or it can some particular

regions where

J-P. Ponthot, T. Belytschko

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

29

this high nodal concentration can be very fruitful, such as around a crack tip in fracture mechanics problem. In such a case, the referential displacement of some nodes around the crack tip can be coupled to the crack tip displacement while other nodes will remain Lagrangian. This technique will be discussed in more detail during the examples. Once the new nodal pattern has been determined, the solution variables have to be remapped with an advection algorithm. Clearly, from (53) to (56), it is seen that new values of the velocity and displacement, at the new node location, are required so that the algorithm can carry over to the next time-step. If the material modeled were path-d.ependent, stresses and stress-associated quantities (i.e. yield stress, back-stresses . .) should also be remapped in the general case where the cells associated with the spatial integration (see Section 3.2) have a relative velocity with respect to the material (see e.g. the wave propagation problem in the numerical examples). However, in some particular cases, it could happen that the spatial integration cells (and therefore the Gauss points associated with those cells) have no relative movement with respect to the material. This is, for example, the case in the two crack-propagation problems (see numerical examples) where only a few nodes located in the neighborhood of the crack tip are declared ALE nodes. All other nodes are located at the vertexes of the integration ~11s and are Lagrangian. Consequently, in that case, all Gauss points are Lagrangian and therefore they require no advection of their associated quantities. In such a case, the extension of the present formulation to accommodate nonlinear path dependent materials (elastic-plastic or elastic-viscoplastic) is rather straightforward and depends only on the availability of a Lagrangian algorithm to integrate those constitutive equations in time. In the ALE/solid finite element community, lots of efforts have been devoted to designing such transport algorithms. These sc:hemes generally rely on the experience gained in computational fluid dynamics. Liu et al. [62-641 extended the SUPG algorithm from Hughes and others [21,44] for Gauss point centered variables. Benson [20,18] generalized Godunov/Van Leer MUSCL (Monotone Upwind Schemes for Conservation Laws) algorithms for both element centered and node centered algorithms. The latter requires the construction of a staggered mesh [20] which is only possible for a structured mesh. H&ink, Akkerman and others [40,39,2,1] presented two algorithms. One is based on an averaging procedure used in post-processing of finite element calculations. The other evaluates the fluxes through the element sides and stems from a finite difference method for compressible fluid dynamics. Both of these schemes are in fact an adaptation of the Lax-Wendroff finite difference scheme. Huerta et al. [37,38] also presented two schemes: one based on a Godunov scheme, the other also relying on a Lax-Wendroff scheme. They also discussed the generalization of the convective algorithms to higher-order Gauss quadrature schemes. Actually, with the exception of the SUPG scheme, all these convective algorithms are built while supposing one Gauss point per element (this is due to the fact that those schemes were generally developed in a finite volume/finite difference community) and their extension to higher-order quadrature is by no way a trivial task as shown in [38]. On the other hand, ALE based fluid-structure interaction codes generally use donor cell upwinding to deal with the convection in the fluid (see e.g. [7,49] or [24-261). Actually, in most schemes, some amount of upwinding is always introduced, by one means or another, in the algorithm. Our problem is now to determine the new value, i.e. the value after convection, of a quantity f (f being any velocity, displacement, stress or stress-related quantity component) on the new configuration. Using (31) with f = 0, since this terrn was accounted for in the Lagrangian phase, one gets

o=f+
(57)

so that we find for a time increment: “f,.

=

j-r+A’f:dt =- l’+A’ (v - ti)*vfdt ,

This equation is valid for both FEM and EFG. However, the FEM solution is tedious since the gradient of nodal values (velocity, displacement) is generally unknown. In fact, since these fields generally only have Co continuity, their gradient is not defined on inter-element boundaries. The same situation prevails for stress and stress-related quantities which are interpolated with C ’ functions whose gradient is not defined. The situation is completely different in EFG because the gradients of the unknown field can be evaluated easily anywhere. Therefore, considering the integrand constant in (58), AflX can be approximated by

30

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 15.2 (1998) 19-46

Aflx=Ar(tj-v).I-$f=d.V

(59)

where d = (ti - u) At is the total displacement during the convective step. Introducing parameter ~9such that 6 e [0, I], the gradient off can in turn be easily evaluated anywhere along the d vector; 19= 0 corresponding to the origin (i.e. the location obtained just after the Lagrangian phase) and 0 = 1 corresponding to the final location (i.e. the location after the convection phase). The final form for any quantity f is thus given by

f

=fLAG + Afl, =fLAG + A@c - ~1.

K=

I,...,n

with XN= XLAG+ 8 At(fi - u)

(61)

where fLAG and xLAGare the quantity f and the position x resulting from the Lagrangian substep. However, due to the fact that the Kronecker delta condition (13) is not satisfied at each node, a distinction must be made between the nodal valuef(x,) and the nodal parameterA at point x,. As the time integration algorithm deals with nodal parameters, relation (60) which results in the nodal value for f have to be modified. Therefore, for node associated quantities, we have adopted the following expression:

(6% i.e. during the transport phase, the nodal parameters are modified according to the gradient of the field quantity. As will be shown in the numerical examples, this approximation seems quite effective. Beyond the fact that, in an EFG formulation (62) is very easy to evaluate for any value of 0, it has the following essential properties: it is consistent, in the sense that for a Lagrangian node, all quantities remain unaffected since, in such a case u = 6 (typical remeshing algorithms for FEM are usually not consistent); it trivially satisfies the DeBar consistency condition (see Benson [20]) which states that if a body has a uniform velocity and a variable mass density before advection, then the body must have the same uniform velocity after advection; it can be applied to nodal quantities no matter how many Gauss points are used in the cell integration process; it does not require the creation of a staggered data structure (i.e. a staggered mesh in a finite element context); it can be applied in 1, 2 and 3 dimensions; in contrast to FEM, where transport is assumed to occur oniy through the cell edges, thus ignoring the coupling between diagonal cells (in a structured 2-D mesh, each element has 8 neighbors but shares its edges-through which a flux is computed-with only 4 of them), EFG advection is carried out into the whole neighbourhood of the node, thus avoiding the loss of accuracy resulting from the lack of (element) comer advection; moreover, as will be shown in the numerical examples, it allows in a natural fashion the introduction of some amount of upwinding, simply by modifying the value of the parameter 19.This is a very important property since, for a long time, upwinding has been considered essential for the good behavior of an advective algorithm. REMARK 3. An alternative way to derive Eq. (57) for the velocity is by using the weak form of the momentum conservation equation while supposing no contribution from internal and external forces (i.e. retaining the first two terms in (37)), so that

Ip

a~,[4 + c,u,,jI

dfi = 0

(63)

R

Choosing 6u to be the Dirac delta function, i.e. a mathematical requirement that the residual of the weak form vanishes at each collocation point, restores relation (57).

J-P. Ponthot, T. Belytschko

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

31

REMARK 4. As shown by relation (63), only the new velocity field can be obtained, after convection, from a conservation equation. The equations for displacement, stresses and stress-related quantities are obtained by analogy with this relation and do not result from a conservation equation. 3.4. ALE algorithm 3.4.1. Pre-processing phase l Define ALE and Lagrangian nodes and their initial position and velocity. l Compute shape functions and their derivatives: - Loop over cells of the domain * Loop over quadrature points xo in cell C * 1. if quadrature point is outside the physical domain, go to 4; 2. check all nodes in cell and surrounding cells to determine k nodes x,, I = 1 to k such that xo is in their domain of influence; 3. for each of the k neighbor nodes, compute @,(x0) and its spatial derivatives; 4. continue; * End quadrature point loop - End cell loop 3.4.2. l

Computational phase Loop over time steps - compute stress u at each quadrature point according - compute internal forces ft”’ = J, BTU da - compute external forces f:“’ - compute acceleration a; = (Mb)-‘(fy”’ -f;“‘) - update velocity, displacement and position

to constitutive

law

,1+112 u,

:=u, n+l

UI LAG XI

l

=

n-“2

u’,’

+

=$’

+

Ata;’

i\tU;+‘12

++‘+“2

- update crack position (if any) - move ALE nodes (if any) llCI x, =.x, LAG+ At(tir - u,) - if crack and/‘or ALE nodes are moving: * update shape functions and derivatives * transfer required variables End loop over time steps

4. Numerical

where required

examples

4.1. One-dimensional

wave propagation

problem

The first numerical example is a one-dimensional wave propagation problem which was introduced by Liu et al. [62] in an ALE/finite element context. This benchmark has an analytic solution and it provides a severe test of the update procedure. Therefore, it has been extensively treated in the FEM literature (see e.g. [62,63,37,38] or [18,20,40] for modified but similar versions). The problem statement, given in Fig. 3 represents a 1-D elastic rod. The spatial discretization is arranged so that no reflected wave will occur during the time interval under consideration thus modeling an infinite rod. Constant density and isothermal conditions are assumed. The rod is modeled by 401 EFG nodes separated by & = 0.1 (all dimensions are in meters). Node 1 is located at the origin of the spatial coordinates and, between

J.-P. Ponthot, T. Belytschko

32

PO)

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

2

1

400

3 l

401

,,-

0.1 W)

40 +

100

I -

time (set)

Fig. 3. Problem description

for 1-D wave propagation.

each pair of EFG nodes is a quadrature cell. One point integration is used in each cell. The equations of motion have been integrated in time with a central difference method. The material properties are the following: E = 10 000 N/m*, v = 0.0, p = 1 kg/m’. These properties are = 100 m/s and the critical time-step for such that the acoustic wave speed in the material is given by c,, =mp FEM is thus 0.001 s, i.e. the time necessary for a pressure wave to cross a cell. The propagation of a square compression stress wave (- 100 N/m’ in amplitude and 0.0045 s in width by the ALE/EFG method) is simulated in three stages: (1) t E [0,0.045] The grid is fixed. A square wave is generated at the origin of the spatial coordinates; (2) t E [0.045,0.08] The grid is fixed and the wave travels along the bar; (3) t E [0.08,0.16] The grid is moved uniformly to the left-hand side with a constant speed ti (= -25 m/s). The resulting stress is then reported at t = 0.16 as a function of spatial coordinates. REMARK

5. The more classical Lagrangian results are simply obtained by letting ti = u in the third step. In fact, in the present case, material displacements are so small that Lagrangian and Eulerian formulations are equivalent to first order; however, the nodal motion has been designed to be large.

REMARK

6. Since the purpose of this simulation is to provide a severe test to the formulation the grid velocity is taken equal to one fourth of the speed of sound.

and the algorithm,

REMARK

7. In [62,63,37], alternative results with a positive mesh velocity are also exhibited. This second case is much less severe than the present one since the relative convective velocity, i.e. ICI = Iu - G 1= 100 - 25, is far less important than the convective velocity considered here which is given by ICI= 100 + 25 = 125. For the same reason, elasto-plastic results are not reproduced here since the plastic effect is to introduce some high-frequency damping which damps out the oscillations (see [62,63]). REMARK

8. For the present problem, two alternative implementations were tested. The first one was a total formulation for the stress-strain relation. This formulation requires the transfer of both the velocity and displacement fields. The second one used an incremental formulation for the stress evaluation and thus required velocity and stress tensor to be convected during the Eulerian phase. Both formulations gave identical results.

Fig. 4 displays the analytical solution as well as the numerical results obtained while using EFG, a linear basis and a Lagrangian description with 320 equal time-steps of 0.0005 s. In FEM, this corresponds to one half of the critical time-step. These results bear a strong resemblance, and in fact they are identical, to the FEM results obtained with linear elements. This simply results from the choice of the domain of influence of the nodes. In fact, the radius of support r, of the nodes is here spatially constant and it can be more easily expressed in terms of the parameter d,, which is defined by rl =

&A,,,

where h,,,, is a characteristic

length of the integration

cells, i.e. in

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. I52 (1998) 19-46

Lagrangian

33

results

25 0 -25 -50 brr -7.5

numerical

-

m1OO -125 -150 spatial

Fig. 4. I-D wave propagation:

.

1-D:

coordinate

m = 2; d,,, =

I; At = 0.0005.

hce,, = LIX

(65) (66)

where Ax and Ay are the cell dimensions. Use of d,, = 1, as is the case here guarantees that each Gauss point has two and only two neighbors in one dimension. Therefore, as already shown by Liu et al. [%I, it can be demonstrated (see Appendix A for an alternative demonstration) that in such a 1-D case, for a linear basis, i.e. m = 2 or p’(x) = [ 1, x], the EFG shape functions are independent of the weight function and are identical to the FEM shape functions. Increasing d,,, (Fig. 5), or the degree of the approximation (Fig. 6) still leads to oscillatory results. However, these oscillations are stable and are completely damped in the neighborhood of x = 0. ALE results for a linear interpolation (m = 2) and d,,, = 1 are shown in Figs. 7 and 8 (all curves are identical) for different values of the parameter 0 from relation (60)-(61). The effect of modifying 0 is striking. For 8 = 0,

Lagrangian

Lagrangian

results

50

results

!

25 0 -25 0..

-50 numerical

-75

-

-100

0

2

4

6

8 spatial

10 12 coordinate

14

16

18

0

20

2

4

6

6 spatial

10 1’2 coordinate

14

16

1.3

20

8 spatial

10 12 coordinate

14

16

18

20

= 2; d,,$= 2; At = 0.0005 s.

Fig. 5. I-D wave propagation:

112

Fig. 6. 1-D wave propagation:

m = 3; d,,, = 2; At = 0.0005 s.

ALE results

75 50 25 0 -25 -50 -75 -100 -125 -150 j

0

2

4

6

8 spatial

’

’

’

10 12 coordinate

14

16

’ 18

20

Fig. 7. 1-D wave propagation: Fig. 8. I-D wave propagation:

0

2

4

6

m = 2; d,>,= I; Ar = 0.0005 s; 0 = 0.0.

m = 2; d,,, = 1; At = 0.0005;

0 = 0.1 and 0 =

1.O.

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

34

i.e. the velocity and displacement gradients are evaluated at their position after the Lagrangian substep, the procedure is clearly unstable and generates oscillations of growing amplitude. This can be understood rather easily, in the present case, since each node J has 3 neighbors in its domain of influence, namely nodes J - 1, J and J + 1. Therefore, the final value of the velocity at a given time step is, according to (60) and (6 1) given by

LAG+ At(v,L*’ - tiJ) $

VJ = VJ

.r=\LAG

with

-au

ax

_!s

I

r=.r)*G

-

and since, according

a@K ax

.I

I

K= J-

,=.,)AC’K

ax

to ((A. 1 l)-see

=.ry =[-1/2AX

l,J,J+

Appendix

1

(68)

“J+;-,“J-I

(70)

A),

1/2hxl

0

one gets

-au ax I

r=.I$AC= [-1/2Ax

1/2Ax]

0

VJ-I VJ

[uJ+I

=

1

i.e. the simple averaging estimation which is known to be unconditionally unstable [21,1]-note the similarity between the EFG relation (70), the finite element relation ((B.4) of Appendix B) and the finite difference relation of Eq. (2.1.7) of [21]. On the other hand, as soon as 8 > 0, the location where the gradient is evaluated is shifted in space and the neighborhood of this point is reduced to 2 neighbors, namely nodes J - 1 and J if ti < v (or J and J + 1 if u^> v) and the EFG estimation of the gradient is given by

-au I ax

a@K

r=x,“=(3x

I

K=J-1,J x=xJHvK

x,”= XJLAC+ 8 At(fiJ

(orK=J,J+l) (72)

- vJ)

where a@Kldx(l,r;lH is given by ((A.lO)-see

au

-I ax

=[_I/&

X=-C:

l/h]

[

Appendix

A) thus resulting

in

“:,’ =vJ-p

(73)

I

which is in fact the upwind evaluation of the velocity gradient (see Appendix B-relation (B.lO)-and relation (2.2.1) of [21]). Thus one can see that simply by choosing 0 to be greater than zero, some amount of upwinding has been introduced into the system of equations. The only overhead work for this inclusion is to recompute the shape functions of the ALE nodes and their derivatives at the location x0, which is straightforward in an EFG code. As can be seen in Fig. 9, results are nearly independent of the time step used. The splitting technique used here to integrate the ALE equations in time allows a critical time step which is, in the present case, 25% higher than the critical time step used in the fully coupled algorithm. The latter is given by [62,25,11]:

Atcrmoupled = fly+ whereas,

ICI = 10::

25 = o.ooo8 s

in the present case, for the Lagrangian

Atcritsplit

=

-=~=O.OOl zp

s

sub step, it is given by the classical

formula

[9]: (75)

This simply results from the fact that during the Lagrangian phase, the convective velocity is equal to zero and that for m = ? and d, = 1, EFG shape functions in 1-D are identical to FEM shape functions.

J.-P. Ponthot, T. Belytschko

I Comput. Methods Appl. Mech.

Engrg. 152 (1998) 19-46

35

ALE results

-25 -50 0.z -75 -100 -

numerical numerical numerical

Cowant Courmt Courant

= 0.5 = 0.7 = 0.9 analytic

4

6

-125 -

0

2

Fig. 9. I-D wave propagation:

8 spatirl

10

12

14

cuordmatr

16

M = 2; d,?,= 1: Ar = 0.0005/7

18

20

and 9 s; 0 = I .O.

For higher d,,, or higher interpolation degree (m > 2) analytical derivation of shape functions is not readily available but the upwind effect is still present as shown in Fig. 10. As also shown on this figure, and contrary to the linear case, the results this time depend on the value adopted for 19. As can be seen by comparing Fig. 9 with Fig. 10, increasing the degree of the interpolation decreases the amount of diffusion. This is confirmed by Fig. 11 which displays the analytical, as well as numerical results obtained for m = 2 and m = 3 at times t = 0.160 (320 time steps) and t = 0.320 (640 time steps). At the final time, the propagating, wave has been convected through 60 cells. 4.2. Mode I crack propagation As already mentioned, ALE formulation is an attractive alternative to the Lagrangian formulation in the sense that it is possible to automatically maintain a dynamically high density of nodes in a limited region evolving with time such as a crack tip. This can be achieved by superposing to a rather coarse Lagrangian discretization a few ALE nodes which move at the same grid velocity ti as the crack tip. One possible way to do this is to generate a few ALE nodes in a star-like configuration (see Fig. 12) around the initial position of the crack tip. Another possibility is to generate a regular node distribution in a rectangular box. As a first example of this technique, we present the following problem of mode I crack propagation. Consider a semi-infinite crack in an infinite body subjected to a pulse of magnitude g, parallel to the crack plane. The body is under plane strain conditions with elastic constants E = 211 GPa, v = 0.30 and p = 7800 kg/m’. The dilatational, shear and Raleigh wavespeeds are denoted as cd, c,, and c,, respectively. At time t = t, , the loading wave reaches the crack plane and at a later time, t = t, the crack starts to propagate at a constant speed u. This problem has been previously analyzed numerically using Lagrangian finite elements by Nakamura et al. [71]. An EFG Lagrangian analysis can also be found in [17,80]. Freund [31] has also investigated an analytical solution for the dynamic stress intensity factor K(t, II). This is given by K(t, u) = k(u)K(t, 0)

(76) ALE results

Fig. 10. 1-D wave: m = 3; d,,, = 2; At = 0.0005; 13= 0.1/0.5 Fig. 11. 1-D wave propagation:

results at time

and 1.O.

t = 0.16 and time t = 0.16; 0 = 1.0.

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

36

0

0

ALE

node

0

Fig. 12. A ‘star-like’

Here, K(t, 0) is the stress intensity &t(

K(t, 0) = 2a,

configuration

for ALE nodes around the crack tip.

factor for a stationary

crack which for this specific case is given by

1 - 2v)ln (77)

1-v

and k(v) is a universal approximated by

function

of crack speed which, for the values of crack speed considered

here, can be

1 - v/c, (78)

k(v) = vz, For the purpose of comparison 1 - Y2 G(t, u) = 7 A(u)[W, where A(v) is also a universal

with numerical

results presented subsequently,

UN* function

the dynamic energy release rate is (79)

of crack speed and is given by

2

A(v) =

’ ad (1 - V)c;D

(80)

Eq. (77) is applicable to crack growth in a bounded body provided reflected waves from the boundaries do not reach the vicinity of the crack tip. To this end, consider a rectangular plate with an edge crack shown in Fig. 13. The dimensions of the plate and crack are chosen such that within the time interval of interest, reflected waves do not approach the crack tip region. The plate is initially unstressed, at rest and completely free. A stress of magnitude a, = 1 is applied as a step function in time to the upper edge in y-direction. Using the central difference method for time integration and a linear basis function in the EFG formulation, the transient response of the plate in plane strain is calculated for 90 time steps with At = 10P5. The dynamic energy release rate G(t, v) is calculated at each time step using the domain integral method [69,70] in a square domain with size = 1.5 m centered at the crack tip. Details of the implementation can be found in [ 171. Solutions were generated with a background Lagrangian grid composed of 100 X 41 (a much coarser grid than the one used in [ 12,17,14,13] where a 200 X 8 1 grid was used) cells with Lagrangian nodes placed at the comers plus 18 ALE nodes in a star-like configuration (2 nodes X 9 angular locations at 0”; 245”; 290”; + 135”; %175”). On each radius of the star, the nodes location was at 0.33331zC,,, and 06666h,,,, where hce,, is the average dimension of a cell according to (66). One point quadrature rule was used in most cells with 6 X 6 Gauss quadrature in regions surrounding the crack. The MLS shape functions were constructed with a linear

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

2.0

31

m

II

4.0 m

5.0m

t

I

10.0

I

m

Fig. 13. Schematic of the discrete model for the infinite plate with a semi-infinite as a step function in time.

crack subjected to a stress wave of magnitude

a,, applied

basis and a Gaussian weight function with a domain of influence based on d,, = 2.0 or d,,, = 2.5. Lagrangian solutions are those reported in [80]. Both Lagrangian and ALE solutions are compared to the analytical solution. Fig. 14 displays the initial grid and crack configuration and Fig. 15 is a magnification of the region surrounding the crack tip. In all these computations and the subsequent, we have used a value of 8 = 1. Three types of solution were generated and compared: (1) a stationary crack, u = 0 (2) a moving crack that begins moving at a constant velocity of u = 0.4c,, immediately when the stress wave hits the crack (3) a crack that remains stationary after the stress wave arrives until time r, and then suddenly begins moving at a constant velocity of u = 0.4c,, The energy release rates for these three cases are shown in Figs. 16 and 17 for d,,, = 2.0 and d,, = 2.5, t, corresponds to respectively. The stress wave arrives at t, = H/c, (2H = 4 m is the height of the plate-time the unit reduced time on the figures), after which the energy release rate begins increasing linearly. The solution is then either along the slope corresponding to u = 0 or u = O.~C,,. For the third case, the crack starts to As can be noted on the figures, the ALE solutions for moving cracks are much propagate at t, = l.:jH/c,. smoother than the L,agrangian ones and ALE results are closer to the analytical solution. The evolution of the crack profile for the second case (u = 0.4c,, immediately when the stress wave arrives) is shown in’Fig. 18; the displacements are magnified by a factor of 10”. The crack tip and the ALE nodes in its neighborhood move from the position x = 5.00 m to x = 5.74 m, moving past 7 nodes. 4.3. Plate with two parallel

edge cracks

Kalthoff and Winkler [45,46] have developed an experimental technique for investigating fracture behavior in high strength maraging steel subjected to in-plane shear loading. The experimental setup consists of a free plate 0

0

0

0.8.0

0

-e-+ l

0

0 1

Fig. 15. Magnified

model for the semi-infinite

view around the crack tip-initial

00

0

0

0

0

0

5

5.1

5.7

ll,llll,,,,,,.

4.Q

Fig. 14. Typical discrete numerical ALE nodes around the crack tip.

0

edge-cracked

configuration.

body. The grid contains

0 = Lagrangian

node;

l

I,

100 X 41 Lagrangian

= Ale node.

nodes plus 18

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

38

2

2 analytic ALE -

1.5

LAG

1.5

-

I

ALE

-

LAG

-

1

6

G 0.5

0.5

0

II

-0.5 0

0.5

I

1.5

2

2.5

-0.5 0

0.5

I

1.5

,

/ analytic ALE LAG

1.5 -

2.5

(a) constant velocily

(a) const~nnf velocity 2

2

%/H

fCd/H

I ti = 0

-

I

2

1.5 -

0

I analytic ALE LAG

I v=lJ

-

0

-0.5 0

0.5

1

1.5

2

-0.5 2.5

kiIH (b) zero velority

followed by constant velocity

0

0.5

I

1.5

2

2.5

WH (b) zwo velocity followed by constant velocity

Fig. 16. Energy release rates for the infinite plate with a semi-infinite by Hr,?,lE. d,,, = 2.0.

crack for the 50 X 21 cell mesh. The energy release rate is normalized

Fig. 17. Energy release rates for the infinite place with a semi-infinite by Ha;lE. d_, = 2.5.

crack for the 50 X 21 cell mesh. The energy release rate is normalized

with two parallel edge notches. The plate, initially at rest and stress-free, is impacted by a projectile as shown in Fig. 19. By changing the acuity of the notch tip and the projectile speed, two regimes of failure were observed for different strain rates. For moderate strain rates, which is modeled here, failure occurs in a brittle fracture mode and the crack propagates at an angle of about 70” with respect to the original notch. At higher strain rates, failure occurs in a shear localization mode due to the generation of a shear band ahead of the notch at a small negative angle (-- 10”). Due to symmetry, only the top half of the plate is modeled in the simulation. The plate dimensions are given in Fig. 19. The body is initially at rest and is subjected to a constant velocity u,, applied on the cracked edge. The impact velocity u0 was chosen by assuming that the section of the plate impacted by the projectile had the same elastic impedance as the projectile itself, so the impact velocity was approximately one-half of the projectile speed [56]. From the experiments by Kalthoff and Winkler [45,46], a representative projectile speed for the brittle fracture mode was 33 m/s, thus u0 = 16.5 m/s. The material properties are for high strength maraging steel (18Ni 1900): E = 190 GPa; v = 0.30 and p = 8000 kg/ m3. To simulate the crack growth in an arbitrary direction, the maximum hoop-stress criterion is used for crack initiation [90]. When the maximum hoop-stress cC,,,treaches a critical value, the crack growth takes place in a direction perpendicular to cCcrit.Once growth is initiated, the crack is restricted to propagate at a constant velocity u = O.l%,, where c, is the shear wave velocity, unless fracture criterion is violated, in which case, crack arrest occurs. The critical hoop-stress is based on the dynamic fracture toughness Kf = 68 MPa&i obtained from mode I experiments. The stress intensity factors KI and K,, are obtained using path-independent contour integrals [69,70]. For mixed mode conditions, the stress intensity factors are decoupled using conservation contour integrals with known auxiliary field [92]. The notch was modeled as a crack. The crack faces were assumed traction-free since the notch faces in the experiments were separated by a finite distance. Thus, in the simulation, the material was allowed to overlap along the crack faces, inducing an initially negative mode I stress-intensity factor.

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

(a) 45 time steps,

h-.,/H

= 1.32

(b) 60 tine

tc.,/H

= 1.76

39

t steps,

i 5omm

(c) 90 time

steps,

ted/H

zoomm

= 2.64

Fig. 18. Evolution of the crack profile for the infinite plate with a semi-infinite crack. Results are for the case 2 (u = 0.4c% at tc,,/H = I) with a 50 X 21 Lagrangian cell mesh +

I8

ALE nodes. The displacements have been magnified by IO”‘.

Fig. 19. Experimental

setup for crack growth in edge-cracked plate.

EFG solutions w’ere obtained with a discretization of 50 X 50 cells with Lagrangian nodes placed at the corners. For ALE simulations, a star-like configuration of 16 nodes, similar to the previous example, was added. The discretization, including ALE nodes and initial crack and J-contour position is shown in Fig. 20. Cell integration was performed with one quadrature point except in the neighborhood of the crack tip where it was automatically increased to 6 X 6 Gauss quadrature as the crack propagated. J-domain integration was computed with 24 X 24 Gauss quadrature. Central difference time integration was used with At = 4.10e7 s until a final time of 80 p,s, i.e. Z!OOtime steps. This is in contrast with previously reported EFG/Lagrangian calculations [17,16,80] which used a 100 X 100 cells grid and a time step of At = 1.10m7 s. The simulation was ended at 80 p,s because shortly after that time, the overlap along the crack faces exceeded the experimental notch width. At that point, the numerical solution began to deviate from the experimental results. To drive the solution further, a contact algorithm should be included. Fig. 21 displays an intermediate (after 100 steps) and the final configuration (200 steps) with the whole crack-path and Fig. 22 is a magnified view around the crack tip. In this figure, one clearly sees that some ALE nodes are very close to Lagrangian nodes (this happens all along the simulation). However, this does not seem to induce any instability. Fig. 23 displays the normalized stress-intensity for Lagrangian and ALE simulations. While the crack is stationary, the stress-intensity factors, which are both negative, increase in magnitude similarly to the analytical

J.-P. Ponthor, T. Belytschko I Cotnput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

40

0’

........................................................................................ ..................................................

.I

..

.....................................

-

.....................

.

......

...........................

:::i::::::::::::::::::::::::::::::::::::::::::::::

1...

.........................

1::::::::::::::::::::::::::::::::::::::::::::::::: ....................................................................... 007s

,::::::::::::::::::::::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::::::::::::::~:::::::: :::::::::::::::::::::::::::::::::::::::::::::::::: :::::::::::::::::::::::i::::::::::::::::::::::::::: _ ................................................. :::::::::::::::::::::::::::::::::::::::::::::::::: ::::::::::::::::;::::::::::::::::::::::::::~::::::

0 c+.:::::::::::::::::::::::::::::::::::::::i:::::::::

i::::::::::::::::::::::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::::::::::::::::::::::::

1::::::::::::::::::::::::::::::::::::::::::::::::: ::::::::i::::::::::::::::::::::::::::::::::::::::: :::::::::::::::::::::::::::::::::::::::::::::::::: _ ........... .* ................ , ............... :::::::::::::::::::::::::::::::::::::::::::::::::: .................................................. M& ,.. ... ‘..... .$ .... . ..... ........

a

0’

. ...

6’~,;1‘o:

oo.,..q.....

‘Ov6~Y’nL,us

drr

IO”

‘ilYk

.................................................. ................................................. :::::::::::::::::::::::::::::::::I::::::~:::~::::: - .................................................

:::::::!:::!!:!:!::::::::::::::::::::::::~::~:! :::::::::::::::::::::::::::::::::::::::::::::::::: . ................................................., , 0015

::::;:::::::::::::::::::::::::::::::::::::::~:: .. .......................................................

-.

.

.....

, .....

....................

......

................................................ .................................................. ..................................................

:::::::::::::::::::::::::::::::::::::::::::::::::: -

..*.

.

......................................

......

;::::::i::::::::::::::::::::::::::::::::::::::::::

1::::::::::::::::::::::::::::::::::::::::::::::::: :::::::::::::::::::::::::::::::::::::::::::::::::~

_

I.. ............................................... .................................................. ................................................. ............. *. ......

. ............................

iiii:i:ii:i::i:iiii:I’ii:lili:i:iI:i::i:~::~:~~~~:

00

0025

/I>,

ODI

ronnguratm alter

c 0015

01

200 stepa

Fig. 20. Typical grid for modeling crack growth in an edge-cracked plate. The discretization contains Sl X 51 Lagrangian nodes + 16 ALE nodes. The initial crack and J-domain are also shown. Fig. 21. Numerical

Fig. 22. Magnification Fig. 23. Stress intensity

crack path in edge-cracked

around crack tip for the edge-cracked

plate problem.

factors during crack growth in an edge-cracked

plate.

0 = Lagrangian

plate obtained

node;

l

= ALE node.

with a 51 X 51 node EFG grid.

solution from Lee and Freund [.%I. Once the crack begins propagating, a different situation occurs: the crack growth is now predominantly mode I. In that figure, one can see that there is a good agreement between both formulations. However, the average slope of the propagating crack is 78” foi the Lagrangian formulation (100 X 100 discretization) and 73” for the ALE formulation (50 X 50 discretization), a value which is closer to the experimental value of 70”.

5. Conclusions A new formulation that combines both advantages of the EFG (easy representation of dynamic crack propagation) and ALE (possibility to automatically keep an evolving high density of nodes where it is required) has been presented. Due to the adopted operator-split methodology, the new formulation is hardly more complicated than the Lagrangian/EFG formulation and can easily be implemented in a pre-existent EFG or

J.-P. ponthot,

T. Belytschko

41

/ Comput. Methods Appl. Mech. Engrg. 15.2 (1998) 19-46

meshless code. A key feature of the new formulation is that some amount of upwinding can, in a natural way, be introduced in the transfer procedure. The method has been applied successfully to wave propagation problems as well as dynamic crack propagation problems. Moreover, if a (rather) coarse Lagrangian background cell is used for the spatial integration purpose, the extension of the formulation to elasto-plasticity and elasto-viscoplasticity only requires a Lagr-angian algorithm to integrate those constitutive laws in time.

Appendix

A. Shape function calculation

for some typical EFG 1-D calculations

The goal of this appendix is to evaluate the 1-D EFG shape function the neighborhood of some point y. A.I.

Case A: 2 nodes in the neighborhood

resulting

from some typical situations

in

of y

This situation is depicted in Fig. A.1 where the weight function W(X) is schematized. Using a linear basis, i.e. m = 2 or pT&) = [I, xl, the A(x) and B(x) of Eqs. (6) and (7) are given by the following matrix form:

(A.11 (A.21 where w, = W(X-x0:1

(A.3)

w*=w(x-X”_AX)

(A.4)

Since P and Ware square invertible

matrices (given the weight function

is strictly positive), it then follows that

-xo

(A.9

1 I so that the shape functions @T=&]l

.X1 [

do not depend on the weight function

x0 +Ax _1

-x0 1 ]=&,x,+iLr-x

and are simply given by (12), i.e.

x-x,J

(‘4.6)

or, by x -

5=

X”

(A.7)

Ax

one obtains @‘=[l-[

[]

(A.8)

Ax Fig. A.l.

Weight function

w(x) such that there are 2 nodes in the neighborhood

of y.

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

42

.

. Y

AX Fig. A.2. Weight function

w(x) such that there are 3 nodes in the neighborhood

i.e. the standard shape functions of the FEM. The derivative of the shape functions calculated from (A.8), or more generally from (12), resulting in (@‘)T = (PO)TIAp’B] + PT[A-‘B]’ and, since in the present case [A-‘B]’

of y.

@’ can also be readily

(A.9) = 0, one finds

(w)T=-&,-ll]

(A.lO)

REMARK 9. The present particular results only hold when the number of neighbors (2) is exactly equal to the number of unknown parameters (m = 2), thus resulting in a square and invertible (though ill conditioned (see [SS])) P matrix. In all other cases, P is a rectangular matrix and the simplification that takes place in (86) is impossible since in that case P -’ is meaningless (however, A is always a square matrix). A.2.

Case B: 3 nodes in the neighborhood

of J

This situation is depicted in Fig. A.2. The main difference from the previous situation is that the radius of the support at point 4’ is now larger and that 3 nodes lie in its neighborhood. Using a linear basis, i.e. m = 2 or pT(x) = [ 1, x]. it can be found that 1 (@‘(X,))T=z]-l

Appendix

B. 1-D averaging

0

(A.1 1)

11

for FEM

Consider the 1-D FEM mesh with nodal velocity u and stress u (1 Gauss point per element) as depicted in Fig. B.1. On each element, the shape functions are simply [ 1 - 5 [I. All elements are identical and have a length of Ax. From the discrete u field, we can obtain a continuous (T” field by a simple averaging:

(B.1) So that on the Jth element,

the stress-gradient

is given by

03.2) where the last equality has been obtained by applying (B.l) twice. From (B.2), it is evident that, according this simple averaging technique, stress gradient on the Jth element does not depend on a,. (T

(JJ

OJ-I

1+1 .

.

J-l

J-2

J

Jcl

J+2

AX Fig. B.I.

1-D FEM mesh with node and Gauss points location.

to

J.-P. Ponthot, T. Belytschko I Cotnput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

The situation element J:

is essentially

the same for the velocity

(Vu), = (VCqJu, = vJ+L where linear shape functions (RI),=

vu,-,-=+vv,

gradient evaluated

at the nodes since one can write on

uJ

(B.3)

have been used. At node J, the velocity

gradient

can thus be estimated

%+I - 'J-l

:!

- a@/+,

~+,(1/2-a)uJ+, where cy is a user-defined 1/2ZaSO -1/2ScuSO

ifc=u-620

gradient,

+ 2q

- (l/2

+ cu)a;~,]

(B.5)

+ 2q

- (l/2

+ a)u,_,]

(B.6)

parameter

ifc=u-tic0

by (B.4)

2A.x

Once again, velocity gradient at node J does not depend on the velocity at that node. This situation can be remedied by introducing some amount of upwinding. The stress and velocity respectively on element J and at node J, are then given by 2-l [(l/2 ax - Ax

43

that has to satisfy the following

constraint (B.7) (B.8)

if (Y= 0, there is nlo upwind and relations (B.2) and (B.4) are recovered. If 1~~1 > 0, some upwind has been introduced. If lcrl = l/2, the algorithm is often termed full upwind or full donor. In this case, for c C 0, the full-donor upwind evaluation of the gradients is given by

aa;__!_(a; - UJ-,I ax - Ax

2 =& (v, .-

(B.9)

VJ_ , )

Acknowledgments Dr. P. Krysl deserves special thanks for all his valuable comments and fruitful discussions regarding many aspects and subtleties of EFG. The support of the Office of Naval Research is gratefully acknowledged.

References HI R. Akkerman, Euler-Lagrange M

131 141 151 161 171 181 191 r101

simulations of non-isothermal viscoelastic flows. Ph.D. Thesis, University of Twente, Enschede, The Netherlands, 1993. R. Akkerman, J. Huttink and P.N. van der Helm, Finite elements and volumes in a Euler Lagrange formulation-artificial dissipation versus limited flux schemes, in: E. Ofiate and D.R.J. Owen, eds., hit. Conf. on Computational Plasticity, Barcelona, Spain (1995) 2307-23 18. J. Argyris, J.-St. Doltsinis, H. Fisher and H. Wurstenberg, TA rr ANTA PEI. Comput. Methods Appl. Mech. Engrg. 51 (1985) 289-362. F.l?T. Baaijens, An U-ALE formulation for 3-D unsteady viscoelastic flow. bit. J. Numer. Methods Engrg. 36 (1993) 1115-1143. I. BabuSka and J.M. Melenk, The partition of unity finite element method, Technical Report BN-1185, Inst. for Phys. SC. and Tech., University of Maryland, Maryland, 1995. I. Babugka and J.M. Melenk, The partition of unity finite element method: Basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289-315. T. Belytschko, D.P. Flanagan and J.M. Kennedy, Finite element methods with user-controlled meshes for fluid-structure interaction, Comput. Methods ,4ppl. Mech. Engrg. 33 (1982) 669-688. T. Belytschko, L. Gu and Y.Y. Lu, Fracture and crack growth by element-free Galerkin methods, Modell. Simul. Mater. Sci. Engrg. 2 (1994) 519-534. T. Belytschko and T.J.R. Hughes, eds., Comput. Methods for Transient Anal. (North Holland, 1983). T. Belytschko and J.M. Kennedy, Computer models for subassembly simulation, Nucl. Engrg. Des. 49( 1) (1978) 17-38.

J.-P. Ponthot, T. Belytschko I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

44

[l 11T. Belytschko, Vessel Technol.

J.M. Kennedy and D.F. Schoeberle, 102 (1980) 62-69.

Quasi-Eulerian

finite element formulation

for fluid-structure

interaction,

J. Pressure

]I21 T. Belytschko, Y. Krongauz, M. Fleming, D. Organ and W.KK. Liu, Smoothing and accelerated computations in the element-free Galerkin method, Comput. Methods Appl. Mech. Engrg. 139 (1996) l-47. 1131 T. Belytschko, Y.Y. LU and L. Gu, Crack propagation by element-free Gale&in methods, Engrg. Fracture Mech. 51(2) (1995) 295-315. [I41T. Belytschko, Y.Y. Lu, L. GU and M. Tabbara, Element-free Galerkin methods for static and dynamic fracture, Int. J. Solids Struct. 32( 17-18) (1995) 2547-2570. 1151 T. Belytschko, [I61 T. Belytschko, 186-195.

Y.Y. Lu and L. Gu, Element-free Galerkin methods, Int. J. Numer. Methods Engrg. 37 (1994) 229-256. D. Organ and Y. Krongauz, A coupled finite element-element-free Galerkin method, Comput. Mech.

17 (1995)

1171 T. Belytschko 923-938.

and M. Tabbara,

39 (1996)

Dynamic

fracture

using element-free

Galerkin

methods,

Int. J. Numer. Methods

Engrg.

[18] D.J. Benson, An efficient, accurate, simple ALE method for nonlinear finite element programs, Comput. Methods Appl. Mech. Engrg. (1989) 305-350. (191 D.J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Comput. Methods Appl. Mech. Engrg. 99 (1992) 235-394. [20] D.J. Benson, Momentum advection on a staggered mesh, J. Comput. Phys. 100 (1992) 143-162. 1211 A.N. Brooks and T.J.R. Hughes, Streamline UpwindlPetrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259. [22] J.P. Cescutti, E. Wey and J.L. Chenot, Finite element calculation of hot forging with continuous remeshing, in: J.L. Chenot and E. Otiate, eds., Modelling of Metal Forming Processes, EUROMECH-233 (Sophia-Antipolis, France, 1988) 207-216. [23] A. Chorin, T.J.R. Hughes, M.F. McCracken and J.E. Marsden, Product formulas and numerical algorithms, Comm Pure Appl. Math. 31 (1978) 205-256. [24] J. Donea, Finite element analysis of transient dynamic fluid-structure interaction, in: J. Don&a, ed., Advanced Structural Dynamics (Advanced Science Publishers, 1978) 255-290. [25] J. Donea, Arbitrary Lagrangian-Eulerian finite element methods, in: T. Belytschko and T.J.R. Hughes, eds., Computational Methods for Transient Analysis, Chap. 10 (North Holland, 1983) 474-516. [26] J. Donea, S. Giuliani and J.P. Halleux, An Arbitrary Lagrangian-Eulerian FEM for transient dynamic fluid-structure interactions, Comput. Methods Appl. Mech. Engrg. 33 (1982) 689-723. [27] C.A. Duarte, A review of some meshless methods to solve partial differential equations, Technical Report 95-06, Texas Institute for Computational and Applied Mathematics, Austin, 1995. [28] C.A. Duarte and J.T. Oden, Hp clouds-A meshless method to solve boundary-value problems, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237-262. [29] M. Fleming, Y.A. Chu, B. Moran and T. Belytschko, Enriched element-free Gale&in methods for singular fields, Int. J. Numer. Methods Engrg, submitted for publication. [30] R.M. Frank and R.B. Lazarus, Mixed Eulerian-Lagrangian method, in: B. Alder et al., eds., Methods in Computational Physics, Fundamental Methods in Hydrodynamics, Vol. 3 (Academic Press, New York, 1964). [31] L.B. Freund, Dynamic Fracture Mechanics (Cambridge University Press, 1990). [32] S. Ghosh and N. Kikuchi, An Arbitrary Lagrangian-Eulerian FEM for large deformation analysis for elastic-viscoplastic solids, Comput. Methods Appl. Mech. Engrg. 86 (1991) 127-188. [33] R.A. Gingold and J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical Sot. 181 (1977) 375-389. [34] S. Giuliani, An algorithm for continuous rezoning of the hydrodynamic grid in Arbitrary Lagrangian-Eulerian

stars, Mon. Not. R. Astr. computer

codes, Nucl.

Engrg. Des. 72(2) (1982) 205-212. [35] R.B. Haber and B.J. Hariandja, Computational strategies for nonlinear and fracture mechanics problems. An Eulerian-Lagrangian finite element approach to large deformation frictional contact, Comput. Struct. 20(1-3) (1985) 193-201. [36] D. Hegen, Element-free Galerkin methods in combination with finite element approaches, Comput. Methods Appl. Mech. Engrg. 135 (1996) 143-166. [37] A. Huerta and F. Casadei, New ALE applications in nonlinear fast transient solid dynamics, Engrg. Comput. 11 (1994) 317-345. [38] A. Huerta, F. Casadei and J. Donea, ALE stress update in transient plasticity problems, in: E. Oilate and D.R.J. Owen, eds., International Conference on Computational Plasticity, Barcelona, Spain ( 1995) 1865- 1876. [39] J. H&ink and P.N. van der Helm, On Euler-Lagrange finite element formulation in forming and fluid problems, in: J.L. Chenot, R.D. Wood and 0. Zienkiewicz, eds, Numerical Methods in Industrial Forming Processes (Sophia-Antipolis, France, 1992) 45-54. [40] J. H&ink and J. van der Lugt, Progress in mixed Eulerian-Lagrangian finite element simulation of forming processes, Int. J. Numer. Methods Engrg. 30 (1990) 1441-1457. [41] T.J.R. Hughes, The Finite Element Method (Prentice-Hall, Englewood Cliffs, New Jersey, 1987). [42] T.J.R. Hughes, Numerical implementation of constitutive models: Rate-independent deviatoric plasticity, in: S. Nemat-Nasser, R.J. Asaro and G.A. Hegemier, eds., Theoretical Foundation for Large-Scale Computations of Non-linear Material Behavior (M. Nijhoff, 1983). [43] T.J.R. Hughes, W.K. Liu and T.K. Zimmerman, Lagrangian Eulerian finite element formulation Comput. Methods Appl. Mech. Engrg. 29 (1981) 329-349. [44] T.J.R. Hughes, M. Mallet and A. Mizukami, A new finite element formulation for computational Comput. Methods Appl. Mech. Engrg. 54 (1986) 341-355.

for incompressible fluid dynamics:

viscous

flows,

II: beyond SUPG,

J.-P. Ponthot, T. Belytschko

I Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

45

[45] J.F. Kalthoff, Transition in the failure behavior of dynamically shear loaded cracks, Appl. Mech. Rev. 43(5) (1990) S247-S250. [46] J.F. Kalthoff and S. Winkler, Failure mode transition at high rates of shear loading, in: C.Y. Chiem, H.D. Kunze and L.W. Meyer, eds., Impact Loading and Dynamic Behaviour of Materials, Vol. I ( 1987) 185-195. [47] E.J. Kansa, Multiquadrics-a scattered data approximation scheme with application to computational fluid dynamics: 1. Surface approximations and partial derivative estimates, Comput. Math. Applic. 19 (1990) 127-145. [48] E.J. Kansa, Multiquadrics-a scattered data approximation scheme with application to computational fluid dynamics: II. Solutions to parabolic, hyperbolic and elliptic partial differential equations, Comput. Math. Applic. I9 ( 1990) 147-161. [49] J.M. Kennedy and T. Belytschko, Theory and application of a finite element method for Arbitrary Lagrangian-Eulerian fluids and structures, Nucl. Engrg. Des. 68 (1981) 129-146. [50] S.W. Key, Current capabilities for simulating the extreme distortion of thin structures subjected to severe impacts, in: A.K. Noor and H.D. Carden, eds., Computational Methods for Crashworthiness (NASA Conference Publication 3223, 1993) 165- 184. [51] H.M. Koh and R.B Haber, Elastodynamic formulation of the Eulerian-Lagrangian kinematic description, J. Appl. Mech. 53 ( 1986) 839-845. [52] H.M. Koh, H.S. Lee and R.B. Haber, Dynamic crack propagation analysis using Eulerian-Lagrangian kinematic descriptions, Comput. Mech. (1988) 141-155. [53] Y. Krongauz and 1’. Belytschko, Enforcement of essential boundary conditions in meshless approximations using finite elements, Comput. Methods Appl. Mech. Engrg. I31 (1996) 133-145. [54] P. Krysl and T. Belytschko, Element-free Galerkin method: Convergence of the continuous and discontinuous shape functions, Comput. Methods Appl. Mech. Engrg., accepted for publication. [55] P. Lancaster and K Salkauskas, Curve and Surface Fitting: An Introduction (Academic Press, London, 1986). [56] Y.J. Lee and L.B. Freund, Fracture initiation due to asymmetric impact loading of an edge cracked plate, J. Appl. Mech. 57: ( 1990) 104-111. [57] T. Liszka and J. Orkisa, The finite difference method at arbitrary irregular grids and its applications in applied mechanics, Comput. Struct. I1 (1980) 83-95. [58] W.K. Liu, Y. Chen, S. Jun. J.S. Chen, T. Belytschko, C. Pan, R.A. Uras and C.T. Chang, Overview and applications of the reproducing kernel particle methods, Arch. Comput. Methods Engrg. 3(l) (1996) 3-80. [59] W.K. Liu, S. Jun. S. Li, J. Adee and T. Belytschko, Reproducing kernel particle methods for structural dynamics, Int. J. Numer. Methods Engrg., submitted for publication. [60] W.K. Liu, S. Li and T. Belytschko, Reproducing least square kernel Galerkin method. (i) Methodology and convergence, Comput. Methods Appl. Mech. Engrg., submitted for publication. [6l] W.K. Liu, Finite element procedures for fluid-structure interactions and application to liquid storage tanks, Nucl. Engrg. Des. 65 (I981 ) 221-238. [62] W.K. Liu, T. Belytschko and H. Chang, An Arbitrary Lagrangian-Eulerian finite element method for path dependent materials, Comput. Methods Appl. Mech. Engrg. 58 (1986) 227-245. [63] W.K. Liu, H. Chang, J.S. Chen and T. Belytschko, Arbitrary Lagrangian Eulerian Petrov-Galerkin finite element for nonlinear continua, Comput. Methods Appl. Mech. Engrg., 68 (1988) 259-310. [64] W.K. Liu, J.S. Chen, T. Belytschko and Y.F. Zhang, Adaptive ALE finite elements with particular reference to external work rate on frictional interface, Comput. Methods Appl. Mech. 93 (1991) 189-216. [65] W.K. Liu and D.C. Ma, Computer implementation aspects for fluid-structure interactions problems, Comput. Methods Appl. Mech. Engrg. 31 (1982) 129-148. [66] Y.Y. Lu, T. Belytschko and L. Gu, A new implementation of the element free Galerkin method, Comput. Methods Appl. Mech. Engrg. 113 (1994) 397-414. [67] Y.Y. Lu, T. Belytschko and M. Tabbara, Element-free Galerkin method for wave propagation and dynamic fracture, Comput. Methods Appl. Mech. Engrg. (1995) 131-153. [68] J.J. Monaghan, An introduction to SPH, Comput. Phys. Comm. 48 (1982) 89-96. [69] B. Moran and CF. Shih, Crack tip and associated domain integrals from momentum and energy balance, Engrg. Fract. Mech. 27(6) (1987) 615-641. [70] B. Moran and C.F. Shih, A general treatment of crack tip contour integrals, Int. J. Fract. 35 (1987) 295-310. [71] T. Nakamura, CF. Shih and L.B. Freund, Computational methods based on an energy integral in dynamic fracture, Int. J. Fract. 27 (1985) 229-243. [72] B. Nayroles, G. Touzot and P. Villon, La methode des elements diffus, C.R. Acad. Sci. Paris, 313, SCrie II (1991) 133-138. [73] B. Nayroles, G. Touzot and P. Villon, L’approximation diffuse, C.R. Acad. Sci. Paris, 313, Strie II (1991) 293-296. [74] B. Nayroles, G. Touzot and P. Villon, Nuages de points et approximation diffuse, in: Seminaire d’Analyse Convexe. Universite de Montpellier 2, 1991. Expose No. 16. [75] B. Nayroles, G. Tcuzot and P. Villon, Generalizing the finite element method: Diffuse approximation and diffuse elements, Comput. Mech. IO (1992) 307-318. [76] W.F. Noh, CEL: a time-dependent two-space-dimensional, coupled Eulerian-Lagrangian code, in: B. Alder et al., eds., Methods in Computational Physics, Fundamental Methods in Hydrodynamics, Vol. 3 (Academic Press, New York, 1964) 117-179. [77] T. Nomura, ALE finite element computations of fluid-structure interaction problems, Comput. Methods Appl. Mech. Engrg. I12 (1994) 291-308. [78] T. Nomura and T..I.R. Hughes, An Arbitrary Lagrangian-Eulerian finite element method for interaction of fluid and a rigid body, Comput. Methods Appl. Mech. Engrg. 95 ( 1992) I 15-138. [79] E. Otiate, S. Idelsohn, O.C. Zienkiewicz, R.L. Taylor and C. Sacco, A stabilized finite point method for analysis of fluid mechanics problems, submitted for publication.

46

J.-P. Ponthot, T. Belytschko

[SO] D. Organ, Numerical University, [Xl]

1 Comput. Methods Appl. Mech. Engrg. 152 (1998) 19-46

solutions to dynamic fracture problems using the element-free

Evanston, Illinois,

J.P. Ponthot, A method to reduce cost of mesh deformation in Eulerian-Lagrangian Modelling of Metal Forming Processes, EUROMECH-233

1821 J.P. Ponthot, The use of the Eulerian-Lagrangian

(Sophia-Antipolis,

method, Ph.D. thesis, Northwestern

[83] J.P. Ponthot, The use of the Eulerian-Lagrangian Zienkiewicz,

formulation, in: J.L. Chenot and E. Otiate, eds..

France, 1988) 65-74.

fern with adaptive mesh. Applications to metal forming simulation, in: D.R.J. Owen

and E. Otiate, eds., International Conference on Computational Chenot, R.D. Wood and 0.

Galerkin

1996.

Plasticity, Barcelona, Spain (1992)

2269-2281.

formulation including contact. Applications to forming simulation via fern, in: J.L.

Numerical

Methods

in Industrial

Forming

Processes (Sophia-Antipolis,

France,

1992)

293-300. [84] J.P. Ponthot, A finite element unified treatment of continuum mechanics for solids submitted to large strains, Ph.D. Thesis, Universiti de Liege, Liege, Belgium,

1995, in French.

[85] J.P. Ponthot and M. Hogge, The use of the Eulerian-Lagrangian

FEM in metal forming applications including contact and adaptive

mesh, in: Chandra and J.N. Reddy, eds., Advances in Finite Deformation ASME

Winter Annual Meeting ASME

(AMD-125,

Problems in Material Processing, Atlanta, Ga., Dec. 1991.

1991) 44-64.

[86] S. Qian and J. Weiss, Wavelet and the numerical solution of partial differential equations, J. Comput. Phys. 106 (1993) [87] P.J.G. Schreurs, F.E. Veldpaus and W.A.M.

Brekelmans,

Simulation of forming processes using the Arbitrary

formulation, Comput. Methods Appl. Mech. Engrg. 58 (1986)

155-l 75.

Eulerian-Lagrangian

19-36.

[88] J.C. Simo and T.J.R. Hughes, General return mapping algorithms for rate-independent plasticity, in: C.S. Desai and E. Krempl, eds., Constitutive Laws for Engineering Materials: Theory and Applications (Elsevier, [89] G. Strang, Linear Algebra and its Applications, [90] D. Swenson and A.R. Ingraffea, Comput. Mech. 3 (1988) [9l]

Modeling

2nd edition (Academic

mixed-mode

Amsterdam,

Press, New York,

1987).

1980).

dynamic crack propagation using finite elements: Theory and applications,

381-397.

T. Yamada and F. Kikuchi,

An Arbitrary

Methods Appl. Mech. Engrg. 102 (1993)

Lagrangian-Eulerian

finite element method for incompressible hyperelasticity,

Comput.

149-177.

[92] J.F. Yau, S.S. Wang and H.T. Corten, A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity, J. Appl. Mech. 47 (1980)

3355341.