Comput Mech DOI 10.1007/s00466-017-1463-7

ORIGINAL PAPER

General framework for dynamic large deformation contact problems based on phantom-node X-FEM P. Broumand1 · A. R. Khoei2

Received: 28 February 2017 / Accepted: 3 August 2017 © Springer-Verlag GmbH Germany 2017

Abstract This paper presents a general framework for modeling dynamic large deformation contact-impact problems based on the phantom-node extended finite element method. The large sliding penalty contact formulation is presented based on a master-slave approach which is implemented within the phantom-node X-FEM and an explicit central difference scheme is used to model the inertial effects. The method is compared with conventional contact X-FEM; advantages, limitations and implementational aspects are also addressed. Several numerical examples are presented to show the robustness and accuracy of the proposed method. Keywords Contact · X-FEM · Phantom-node · Large deformation

1 Introduction Contact mechanics is an important branch of mechanics which studies the interaction of solid bodies; it has many application in various fields of science and engineering such as metal forming [1], crash worthiness analysis [2], perforation [3], biology [4], geology [5], fracture mechanics [6], etc. Historically, the finite element method (FEM) has been a common tool to numerically model the contact and impact

B

P. Broumand [email protected] A. R. Khoei [email protected]

1

Department of Civil Engineering, Shahrood University of Technology, Shahrood, Iran

2

Department of Civil Engineering, Center of Excellence in Structures and Earthquake Engineering, Sharif University of Technology, P.O. Box 11365 9313, Tehran, Iran

problems. In fact, modeling the real world problems consists of two steps; first a CAD model is made from the geometry of the prototype and second, suitable meshes are produced based on the CAD model. However, the finite element method requires that the meshes conform to the boundaries, interfaces and existing discontinuities; thus, the conversion of a complicated CAD model into a robust FEM mesh, if not impossible, is very expensive and time consuming. The problem deteriorates when one deals with propagating interfaces like cracks which requires remeshing steps; in addition the FEM is weak in capturing the singularities and the high gradients in the strain and stress fields which are very common in contact and fracture problems. Isogeometric method [7] is an attempt to fill the gap and to use the CAD model directly in the numerical analysis. The method has recently become very popular but still needs more development to become a mature tool for practical uses. In this paper, we wish to propose an alternative framework based on phantom-node X-FEM. The extended finite element method (X-FEM) [8] was first introduced to handle the difficulties with FEM. Since its introduction, the X-FEM has been the subject of several investigations in the field of contact mechanics. For instance, Dolbow et al. [6], Khoei et al. [9], Khoei and Nikbakht [10], Elguedje et al. [11], Liu and Borja [12,13] and Hirmand et al. [14], have implemented frictional contact in small deformation regimes within the 2D X-FEM framework. Geniaut et al. [15] extended the contact formulation to 3D problems. Several attempts have been done to model large sliding contact; for example, Nistore et al. [16] employed a master-slave approach to model static contact problems with large sliding; Liu and Borja [17] used X-FEM to model finite deformation embeded frictional cracks and, Khoei and Mousavi [18] developed a penalty contact method for large sliding problems.

123

Comput Mech

A recently developed version of X-FEM is Phantom-Node method which was introduced by Song et al. [19] to model dynamic shear band and crack propagation. The method has some advantages which makes it attractive for contact problems. First, unlike the X-FEM it is easy to implement in existing finite element codes; second, the method separates the kinematics of the displacement field in the both sides of an interface and thus removes the leaking stress waves from one side to another which are observed in classic X-FEM formulations [20]. The method has been applied to shells with arbitrary cracks by Chau-Dinh et al. [21] and it has become so popular that it has been implemented within commercial finite element toolboxes like ABAQUS. In this paper a general framework for contact and impact problems based on phantom-node method is introduced. For simplicity, the method is proposed for 2D cases but an extension to 3D is almost straightforward. In this method, with help of phantom-Node X-FEM, the geometry which can be obtained from a CAD software is defined independent of the mesh; thus simple structured meshes can be used for very complicated problems. The contact constraint is imposed with a penalty method based on a master-slave approach within the phantom-node formulation. Explicit central difference method is employed with an updated Lagrangian large deformation, large strain approach for temporal and spatial discretization. The paper is organized as follows; in Sect. 2, the generic contact problem is described and the governing equations and the discretization methods are explained. Section 3 is devoted to the proposed general framework for contact problems; several details and obstacles are addressed. Finally, numerical examples are presented in Sect. 4 to verify the accuracy and robustness of the method.

2 The generic contact problem 2.1 Governing equations Without loss of generality, consider a body  bounded externally by surface ∂ with an internal interface , as shown in Fig. 1. The interface divides the body into two separate seg − + and − with the properties that  = +  ments,  and + − =Ø. In addition, some parts of the interface may experience contact which are denoted by c+ and c− corresponding to the either sides of the interface. At time t = 0 the bodies are in the reference configuration 0 C and they are subjected to the initial and boundary conditions. Under the loading and boundary conditions, they deform and occupy the new configuration t C. The motion of a particle is described at the reference position 0 x and the current position t x at time t with the mapping t x = φ(0 x, t). The equation of motion of the current configuration at time t can be expressed as

123

Fig. 1 Contact problem definition

ρ u¨ = ∇.t σ + t b in 

t t

(1)

which is subjected to the boundary and initial conditions as 

0

   x, t = t ud 0 x, t

t

u

t

σ .n = t

t

σ .nc = tc

t

σ .(−n)c = −t tc

on ∂u on ∂t

t

on c−

t

on c+

0

u(0 x, 0) = u0 (0 x)

in 

0

˙ x, 0) = u˙ 0 ( x) u(

in 

0

0

(2)

where tρ is the density, t u is the displacement vector, t σ is the Cauchy stress tensor, t b is the body force vector, t t is the surface traction, n is the normal to the body surface, tc and nc are the contact surface traction and the normal vector respectively. Here (∇.) stands for the divergence operator; the derivatives with respect to the time are also shown with dotted symbols. We show vectors and tensors with boldface fonts and scalars with ordinary fonts. Equation (1) is non-linear in nature; one way to handle it is to use a dynamic explicit formulation which is fully compatible with an incremental Updated Lagrangian approach. In such a case, the equation of motion is considered at time step t and all variables are determined at time step t + t directly; hence, it is not necessary to use a nonlinear solution method like the Newton-Raphson procedure. Since t+t x = 0 x + t+t u and t x = 0 x + t u , we have t+t x = t x + t u and u = t+t u − t u. As Fig. 2 shows, t by applying the virtual work principle to the known configuration at time t and employing the discontinuous version of Gauss–Green Divergence theorem, one leads to the following weak form of the equation of motion

Comput Mech

Fig. 2 Updated Lagrangian approach in dynamic explicit formulation; the displacement variation is applied to the known body configuration at time t and the body configuration at time t + t is found from temporal disceretization

 (δt uT )tρ t u¨ d t V   T t t =− (δt e ) σ d V − δt uT t tc d tS t c−   T t t (δt u ) b d V + (δt uT )t t d tS +

elastic strain E and the Kirchoff stress τ as

t

t

(3)

∂ t t

where, δt u is the virtual displacement imposed at time t, and δ t e = δ(∇ s t u) is the corresponding virtual strain and ∇ s is the symmetric gradient operator. The second term in the right hand side of the equation is the contribution of the contact force; here, δu is the jump of δu across the discontinuity and t tc = t σ .nc− is corresponding contact traction. The above formulation may seem like the small deformation formulation, but here the derivatives are taken with respect to coordinates at time t and all integrals are calculated with respect to the configuration at this time. The Updated Lagrangian formulation presented above, accounts for large deformations, however in many problems the strains are also large. For this reason a hyper-elastic large strain formulation based on reference [22] is implemented. The elastic Hencky strain which is conceptually an extension of the true strain in uniaxial loading into the multi-axial case, is defined as E=

1 1 ln (FFT ) = ln B 2 2

(4)

where F is the elastic deformation gradient, B is the left Cauchy-Green tensor and ln(.) is the tensorial logarithm function. The existence of the Hencky stored energy function results in a linear relationship between the Hencky

ρ =

1 ECE 2

and

τ =ρ

∂ = CE ∂E

(5)

where C = 2μIs + (K − 23 μ)I ⊗ I is the forth order elastic properties tensor, Is is the forth order symmetric projection tensor, I is the second order identity tensor and, μ and K are the shear and bulk modulus. The stress update procedure can be described as following. After obtaining the displacement field t+t u from the finite element solution, since t+t x = 0 x+ t+t u, the deformation t+t t+t gradient is found as t+t F = ∂ ∂ 0 x x = I + ∂ ∂ 0 x u . The Hencky strain is calculated by Eq. (4). In order to compute the tensor logarithmic B, we use the Taylor’s expansion of the function ln(I + X) around the zero tensor as following 1 1 ln(I + X) = I + X − X2 + X3 + · · · 2 3

(6)

where X is a second order tensor. The logarithm of tensor B can be obtained by setting X = B − I. For accuracy, we use the first 200 terms of the series. Equation (5) gives the Kirchhoff stress t+t τ which is related to the Cauchy stress t+t through the relation t+t σ = t+t τJ , where J is the Jacobian of the deformation gradient. 2.2 Discretizing the weak form of the governing equations The weak form (3) needs to be discretized. The original formulation of X-FEM assumes a discontinuous displacement field based on the enrichment of the finite element space by

123

Comput Mech

the partition of unity method [8]. 

u(x) =

I ∈I std



N I (x) d I +

standard part



 I ∈I H



N I (x)H (x)a I

(7)



enriched part

where, H is the Heaviside function defined as H (x) =

−1 x ≤ 0 1

x >0

(8)

Istd and and I H are the set of standard and enriched nodes respectively. d I and a I are the corresponding standard and Heaviside nodal degrees of freedom. For more details on XFEM consult [23]. In dynamic explicit problems, the classic X-FEM discretization suffers from an intrinsic deficiency which can easily be understood with the following wave propagation problem. Similar to reference [24], consider a block with an interface in the middle as shown in Fig. 3. The dimensions are h = 1 cm, w = 2 cm and material properties are those of PMMA: Youngs modulus E = 3.24 GPa, the Poisson ratio ν = 0.35 and density ρ = 1190 kg/m3 . A prescribed velocity function as below, is applied to the top face of the block, V0 t/tr t ≤ tr V (t) = (9) V0 t > tr where tr = 0.1 µs and V0 = 10 m/s Application of the prescribed velocity produces a step wave which propagates through the block. Analytically, it is expected that the wave does not pass through the interface into the other side. However, as Fig. 4 shows, the classic XFEM formulation with mass lumping, transits the wave into the other side of the interface. In other words, in the X-FEM formulation, mass lumping suppresses the correct coupling of the regular and enriched DOFs, thus the two sides of the cut element are kinematically dependent and waves from one side leak to the other side. This is in agreement with the observations of Nistore et al. [20] and de Borst et al. [25] as well. In addition, even if the special lumping technique proposed by Menouillard et al. [24] is used, the problem exists. Although, consistent mass matrix solves this problem, it has two drawbacks: One, using mass lumping is crucial for explicit time stepping; since, otherwise solution of system of equations becomes necessary which increases the computational cost drastically; two, consistent mass matrix may become singular in case of the interfaces which are very close to the nodes. In order to tackle these issues, we employ the phantomXFEM. This method is a variant of the X-FEM that was first proposed by Song et al. [19] to model shear band localization and crack body discontinuity. The method is based on the works of Hansbo and Hansbo [26]. In this method, instead of

123

Fig. 3 A block with an internal interface subjected to prescribed velocity in the top face

enriching the nodes of an element which is cut by an interface with Heaviside functions, two elements are superposed; each element has some real and some phantom nodes. As Fig. 5 shows when an interface completely cuts an element, the displacement field is calculated based on the summation of the displacement field of the two elements; one is u1 (x) which holds for φ(x) < 0 and the other is u2 (x) for φ(x) > 0, u(x) =

 I ∈I1



N I (x)H (−φ(x))u1I +

u1 (x)



 I ∈I2



N I (x)H (φ(x))u2I

u2 (x)

 (10)

where, φ is the signed distance function which separates the two sides of the interface and, 1 and 2 are the sets of the nodes of element 1 and 2 respectively. In order to model discontinuity in the displacement field and stress concentration in crack tip region in case of cracks, different sets of phantom nodes must be introduced. Figure 6 shows different cases that occur in contact and fracture problems. Considering the previous wave propagation problem, it is observed from Fig. 4 that unlike the X-FEM, PhantomNode method performs well and no wave leaks into the other side of the interface. In fact, phantom node method completely separates the two sides of the interface. In case of contact problems this property is crucial. Using the phantom node displacement field (10), the discretized form of the Eq. 3 at time t is obtained as below, t

Mt u¨ = t Fint − t Fext − t G

(11)

since, each element cut by the interface is defined as the superposition of two separate elements e = 1 and e = 2, the element matrices can be expressed as,

Comput Mech

Fig. 4 The σ yy stress value in the block height “y” at a t = 2 µs, b t = 3 µs and c t = 4 µs

Fig. 5 In phantom node method the element cut by an interface is divided into two elements with real and phantom nodes [19]

123

Comput Mech

Fig. 6 Different cases for phantom and real nodes, a an interface cuts the element into two sections, b a crack cuts the element into two section but the crack tip is at the right edge and c a crack cuts two neighboring edges and the crack tip is on the lower edge

 t

M = e

t e Fint t e Fext



t



t

= =

ρNT NH ((−1)e φ(t x)) d tV

(12)

BT t σ H ((−1)e φ(t x)) d tV

(13)

e

Mt+t u = (t)2 (t Fint −t Fext −t G) + M(2(t u) −t−t u)

e

∂t 

NT t t H ((−1)e φ(t x)) d tS

(16)

e

+  t

G=

c−

t

NT t bH ((−1)e φ(t x)) d tS

(14)

t

u¨ =

e

NT t tc d tS

(15)

where, M is the mass matrix, N is the matrix of shape functions, Fint is the vector of internal forces, B is the matrix of the gradients of the shape functions like the ordinary FEM, σ is the vector of Cauchy stress, Fext is the vector of external forces, G is the vector of contact forces, NT  is the jump in shape function matrix and tc is contact traction vector. The equation of motion (11) is a system of second order differential equations that must be integrated in time. In this study, the central difference method is employed, which is an explicit direct time integration method with the second-order accuracy. In comparison to the dynamic implicit method, the dynamic explicit method has several advantages; it is simple, needs a low storage space, does not need to solve the nonlinear system of equations iteratively, and it does not require the stiffness matrix assembly. The dynamic explicit method is the most suitable technique for transient (short time) problems. The method can be used for quasi-static analysis of problems with highly non-linear behavior like contact, that usually have convergence issues. The main disadvantage of the method is that it is conditionally stable. The stability criterion enforces the small time step t particularly for the fine and distorted meshes. Explicit-implicit time stepping methods like the one proposed by Rabczuk and Belytschko [27] are a remedy to this issue. In central difference method, the

123

set of discretized equations in time and space can be written as

t

u˙ =

t+t u

− 2(t u) +t−t

u

(17)

(t)2 −t−t u 2t

t+t u

(18)

The method is not self starting, hence, it needs two conditions for the start of the solution as −t 0

(t)2 0 ¨ ( u) 2 −0 Fext −0 G)

u =0 u − t 0 u˙ +

u¨ = M−1 (0 Fint

(19) (20)

The stability limit tc is approximately equal to the time of an elastic dilatational (pressure) wave to cross the smallest dimension of the model lmin , i.e.

lmin tc = cd

and

cd =

λ + 2μ ρ

(21)

where cd is the speed of dilatational wave and, λ and μ are the Lame’s constants. If an appropriate lumped mass matrix is employed, the central difference method can be used efficiently [28]. In this paper, we assume a simple mass lumping strategy which distributes the mass uniformly among the real and phantom nodes based on the row summation method; this makes the stability limit insensitive to the proximity of the interface to the nodes.

Comput Mech

3 General framework for computational contact problems 3.1 Modeling contact with phantom-node X-FEM Since, the contact is accompanied with large deformations and large sliding, the contact area changes with time and different interfaces may experience contact; thus, the small deformation formulations like, [12] are no longer valid; in order to compute the contact term (see Eq. 15), we need to employ a master-slave scheme. In this method, the contact constraint is applied between a point and a surface, similar to node to surface method in ordinary FEM. The element which contains the point is slave and the element which contains the surface is the master. We assume a criteria for the selection of the master and slave elements: when ever an element is cut with an interface, two elements are produced with the phantom and real node sets as shown in Fig. 6. We take the first element, in which the real nodes have positive level set function values i.e. φ(x) > 0, as the master element and the second element in which the real nodes have negative sign, i.e. φ(x) as the slave element. Considering, Gauss integration method, the Gauss points which are placed on the interface in the slave element are the slave points and the interface line in the master element is the master surface. The interface is defined on the initial configuration and it is updated in every time step. In fact, every contact interface segment (master or slave) is modeled as a line in an element which is cut by the interface (Fig. 7). The interface location is updated by updating the coordinates of the two intersection points of the interface line with the element edges, according to the below formula. t+t

x=

4 

N I (ξ, η)(0 x I +t+t u)

(22)

I =1

Fig. 7 A contact interface is defined on the initial configuration and it is approximated with a line in each element. The two intersection points of each line segment are used to update the coordinates of the interface in subsequent times

Fig. 8 The penalty contact method; a slave integration point G is in contact with a master line. The projection of the integration point on the interface line is the point Gpr

where ξ and η are the local coordinates of the intersection points in the corresponding element. The contact impermeability constraint can be applied with several methods e.g. Lagrange multipliers, Augmented Lagrange multipliers, penalty method, Nitsche method, etc. We assume a penalty approach which is easier to implement and has lower computational cost compared to the other methods. The penalty method assumes that there is a spring between slave point and the master surface (Fig. 8). Any penetration of slave point into the master surface is penalized with a traction force which is proportional with the normal gap. The normal gap is distance between the Gauss point G and corresponding projection point on the master line G Pr ; it is thus computed as

where, ξG and ηG are the parent coordinates of Gauss point G in the slave element; ξG Pr and ηG Pr are the parent coordinates of G Pr in the master element and x I is the nodal coordinates vector of node I . The corresponding normal contact traction on the slave element is thus,

gn = (xG Pr − xG ).nc−

tcn = tcn nc−

(23)

where (.) is the scalar vector product operator and xG = xG Pr =

slave 

N I (ξG , ηG )x I

I master 

and

N I (ξG Pr , ηG Pr )x I

(24)

I

and

tcn = n gn

(25)

123

Comput Mech

where n is the normal penalty coefficient. The contact traction is subjected to the Kuhn-Tucker inequalities as below gn ≥ 0,

tcn

≤ 0,

gn tcn

=0

(26)

Introduction of Eq. 25 into 15 yields,  t

G=

c−

n NT nc− (nc− )T N t x d tS

(27)

for 2D four noded bilinear elements, the jump of the shape function matrix N can be expressed as,  N =

 N1m 0 . . . N4m 0 −N1s 0 . . . −N4s 0 0 N1m . . . 0 N4m 0 −N1s . . . 0 −N4s (28)

where, Nim and Nis are the shape functions of the ith node in master and slave elements. Note that Nim , are calculated in the position of G Pr point and Nis are calculated in the position of G points. Equivalently, the t x vector is defined as below, t

x=

t

x1m t y1m . . . t x4m t y4m t x1s t y1s . . . t x4s t y4s

T

(29)

where, t xim , t yim , t xis and t yis are the x and y coordinates of the ith node of the master and slave elements respectively. In large sliding motions, elements which are originally close to each other and in contact may undergo large sliding and become completely apart and new elements come into contact. For instance, consider Fig. 9, the two elements which are labeled E1 originally coincided; the green element was slave and the blue one was master, but due to motions and deformations, now master element E1 is in contact with E4 and slave element E1 is in contact with E2. A search algorithm is implemented to determine that each slave integration point is in contact with which master interface line at each time.

3.2 General framework for contact The process of converting a physical problem into a discrete mathematic model with help of finite element can be cumbersome. If the geometry of the problem is complex, sophisticated meshing step is required. The issues becomes more complicated in explicit dynamic problems where the size and aspect ratio of the elements has significant impact on the critical time step and thus the computational costs. It is thus very appealing to separate the meshing step from the definition of the geometry. The geometry can be obtained from a CAD software which defines the interfaces that shape the boundaries of the bodies in a problem; the meshing is also preferred to be as simple and structured as possible. The method presented in the previous section has the ability to form a framework which is suitable for this purpose. Consider that we want to model a problem like that of Fig. 10. The bodies are general and may undergo contact during the motion. In the first step, a simple mesh is constructed. The body geometry is then defined with help of interfaces (Fig. 10a). The interfaces may constitute all or part of the boundaries of the bodies. Next, the phantom-node method is applied to each body. Two cases may occur; first, the interface separates the interior of a body from null space; in this case, the elements which are cut by the interface are not doubled but are converted to elements with phantom and real nodes (Fig. 10b). Second, the interface is within a single body, the elements which are cut by the interface are converted to two elements as it was stated in Sect. 2.2 with the corresponding real and phantom nodes and material properties. Finally, all elements which are not in the material domain are excluded from the model problem and the bodies are labeled as master and slave which will be used in the contact calculations. An important remaining issue is how to apply the tractions (Neumann conditions) and Drichlete conditions to the interfaces which are not necessarily matched with the mesh. Tractions on the interfaces can easily be applied with Eq. 14.

Fig. 9 The master slave contact scheme; Gauss integration points in the slave element are the slave points and interface lines in master element are the master surfaces. Due to the large sliding, slave points may experience different masters

123

Comput Mech

Fig. 10 The general framework for contact problems a complicated bodies can be easily modeled with interfaces in simple structured meshes b body meshes are formed with help of phantom and real nodes

Fig. 11 Node relocation; a node in position A which is very close to the interface is relocated to a new position like A’ in a direction normal to the interface

One way to apply the Drichlete boundary conditions weakly, is to use a penalty method. The procedure is the same as contact constraint application, with the difference that no Kuhn-Tucker condition is required. 3.3 Numerical integration In the equation of motion 11, the computational cost is mainly due to the calculation of the internal force term t Fint which must be performed numerically in the Gauss integration points. Reducing the number of Gauss points not only reduces the calculation cost of the integrals but also reduces the number of solutions required for the expensive material constitutive equations. For this purpose, the reduced integration method that is the most popular technique in an explicit

Fig. 12 The geometry and boundary condition of a block with an interface subjected to a downward motion in contact

formulation is employed using only one integration point for a Q4 bilinear elements. The method is also a remedy for the locking issue in lower-order elements of incompressible and nearly incompressible media. However, the reduced integration method results in the singular element stiffness matrix and extra zero energy modes called as the hourglass modes. These zero energy modes do not affect the strain and stress fields but can arouse the displacements significantly. The concentrated forces and restrained boundary conditions are more susceptible to excite these modes. Various techniques have been proposed to control the hourglass phenomenon; among them are the artificial damping and artificial stiffness methods introduced by Flanagan and Belytschko [29] which we used in standard elements.

123

Comput Mech

Song et al. [19] proposed a reduced integration method for phantom-node elements which is very attractive and works very well for traction free interfaces, however our observations showed that this method leads to unrealistic stress and strain fields for contact problems. In fact it was necessary to integrate the finite element equations exactly in these elements. There are two common integration methods in the literature for elements with internal interfaces, namely, the triangulation and rectangulation integration methods. We consider triangulation method in which each side of a cut element and its phantom counterpart is divided into several triangles and integration Gauss points are placed within each triangle.

3.4 Obstacles and limitations Different issues arise when the phantom node X-FEM is accompanied with large deformation contact. We call the most significant one, the hanging node problem. The problem occurs when a contact interface is very close to a mesh node. In this case the master or slave element (depending on the interface location) may have a very small surface area. When, contact tractions which can be very large, are applied to this element, some phantom nodes may undergo unfavorable movements and oscillations; in addition these nodes considerably affect the stress field and contact traction quality. These hanging nodes are very common in problems with

Fig. 13 The σ yy stress contour for a block with an interface at t = 0.1 s a Unstructured coarse mesh, b unstructured fine mesh and c 11 × 11 mesh

Fig. 14 Contours of a u x and b u y displacement for a block with and interface

123

Comput Mech

Fig. 15 Contact pressure profile at time t = 0.1 s for a 11 × 11, 21 × 21 and 31 × 31 structured and two unstructured meshes and b the 11 × 11 mesh with different contact penalty factors

unstructured meshes or in structured meshes with oblique interfaces. The same problem exists in ordinary X-FEM; the usual remedy is to remove the Heaviside enrichment of the nodes which have a small support. This strategy can be translated into phantom node X-FEM as following: the original nodes which do not have big support will not produce phantom nodes. In this manner the produced master and slave elements even without existence of crack tip will have common original nodes and are connected. The problem with this strategy is that in case of large sliding contact, if the master and slave elements are connected, the elements are stretched which apart from its bad view, can lead to large element ratios and low critical time steps. One good remedy for hanging node issue is node relocation. As Fig. 11 shows, in this method, a node, in position A, which is very close to an interface is relocated to a new position like A’, in the direction normal to the interface. The process is needed to be done just before the analysis and the construction of master and slave elements; hence, it does not produce extra computational overhead. Note that in some cases, this simple node relocation method deteriorates the elements’ ratios, hence in general a more complicated algorithm might be necessary to control the elements’ quality. Another remedy which can fully remove the hanging node problem is to artificially increase the internal force in the elements with this problem. For this purpose, one can apply a coefficient α > 1 to the internal force vector (13) as  t e Fint = α BT t σ H ((−1)e φ(t x)) d t V (30) t

e

This will increase the internal force vector value but does not change its direction, thus it increases the stiffness of the

Fig. 16 The geometry and boundary condition of the two body contact problem

deficient element. Note that this artificial force violates the total balance of the system, but since these elements are very small, its impact on the total solution is negligible.

4 Numerical examples 4.1 A block with a horizontal internal interface in contact The first numerical example is presented to verify the accuracy and convergence behavior of the proposed contact method. The problem is a simple 1 m×1 m rectangular block with an interface in the middle. The material has a Young’s

123

Comput Mech

Fig. 17 The σ yy stress contour at time a, b t = 0.1 s, c, d t = 0.15 s and e, f t = 0.2 s; left column ABAQUS, right column this formulation

modulus of E = 10 GPa, Poisson’s ratio of ν = 0.3 and density of ρ = 7800 kg/m3 . The geometry and boundary conditions are shown in Fig. 12. The bottom face is fixed clamped and the top face is fixed in x direction and is sub-

123

jected to a downward displacement of 0.1 m in y direction which is applied by a ramp function in 0.1 s. The loading is quasi-static and hence the results are in steady state condition. The problem is in plane strain state and it is solved

Comput Mech

Fig. 18 von-Mises stress contour at time t = 0.1 s for a 5 × 20 mesh, b 5 × 40 mesh, c 5 × 60 mesh and d 5 × 80 mesh

for five different mesh densities; the first three meshes are structured with 11 × 11, 21 × 21 and 31 × 31 elements, the last two meshes are unstructured and consist of 614 and 2872 degrees of freedom respectively. The contact integration is performed with two Gauss points. A normal penalty factor of 100E is assumed and the time step is chosen to be 10−6 s. Similar problem has been investigate by Liu and Borja [13] for small deformation contact cases. Figure 13 show the σ yy stress contour for the two unstructured meshes and 11 × 11 structured mesh. It is observed that in all cases the stress contours are smooth and almost identical. In addition, the tractions are continuous along the interface as it was expected. Figure 14 shows u x and u y displacement contours at time t = 0.1 s for 11 × 11 mesh. Due to the contact pressure, the block dilatates in the middle. Figure 15a shows the contact pressure profile along the interface for the different meshes. It is evident that the results converge to a unique solution. It is also observed that there are some slight oscillations in the unstructured meshes. These oscillations occur mainly in elements in which the interface is very close to nodes and are prone to hanging node issue. These effects become less significant as the mesh is refined. The effect of penalty contact factor on 11×11 mesh is investigated

in Fig. 15b. Three penalty factors of n = 100E, 1000E and 10000E are considered; it is observed that the solutions are not very sensitive to the penalty factor. 4.2 Two elastic body contact The purpose of this example is to study the performance of the proposed method in the contact of two elastic bodies. As Fig. 16 suggests, the problem consists of two identical curved bodies with internal radius of R1 = 0.25 m and external radius of R2 = 0.49 m; the upper body is pushed for 0.1 m against the lower body and then it is dragged for 0.5 m in the horizontal direction. Both bodies have the same material properties as the previous example and the state of stress is plane strain. The problem is solved for four mesh densities of 5 × 20, 5 × 40, 5 × 60 and 5 × 80 elements respectively. Since the contact area is of greater importance, the refinement is performed only in the tangential direction. The gray lines show the contact interfaces, blue circles are the contact Gauss points and red circles are the projection points of the slave Gauss points on the master interface. A normal penalty factor of 100E is also assumed. As Fig. 16 shows, the loading is displacement controlled and it is applied to the restrains of

123

Comput Mech

Fig. 19 The contact pressure profile on the upper body interface at different time steps a, b 5 × 20 mesh, b, c 5 × 40 mesh, d, e 5 × 60 mesh and f, g 5 × 80 mesh

123

Comput Mech

Fig. 19 continued

the upper body in two phases according to the following functions. u(t) = v(t) =

0 t ≤ 0.1 5(t − 0.1) 0.1 < t ≤ 0.2 t

t ≤ 0.1

0 0.1 < t ≤ 0.2

Figure 17 shows the σ yy stress contours for the 5 × 20 mesh at different time steps. The results are in good agreement with that of ABAQUS-explicit general contact. In the first phase, as the upper body comes down, the contact area grows in size and the stresses increase. In the second phase, when the upper body is dragged against the lower one in the horizontal direction, the contact area moves to the right and the bodies deform anti-symmetrically. Figure 18 compares the von-Mises stress contour at t = 0.1 s for the four mesh densities. It is observed that the results almost converge to a unique solution with mesh refinement. Figure 19 shows the contact pressure profile on the upper body (slave) interface, at different time steps for the four mesh densities. The first column shows the contact pressure during the first phase in which the two bodies are pushed against each other and the second column shows the pressure profile during dragging phase. Considering 5 × 20 mesh, some bumps are observed in the final stages of phase one and two. These bumps occur in elements in which one Gauss point is in contact and the other is not. This effect is almost removed with further refinement of the mesh. As [16] suggests using Newton-Cotes integration method may also alleviates these oscillations. The figure shows that in the first phase of loading, with downward movement of the upper body, almost in all cases, the contact pressure increases smoothly and the

contact area grows in size. In the second phase as the upper body is dragged to right the contact area moves and its size decrease. The results in some instances are less smooth in this phase but are almost identical for different meshes. In order to investigate the effect of penalty factor, the problem with mesh 5 × 20 is solved for three different penalty factors of n = 10E, 100E and 1000E. Figure 20 shows the σ yy stress contour. It is observed that the stress distributions are almost the same; however, for n = 10E the deformed shape shows some penetration which is not favorable, but in the other two cases there is no visible penetration. Figure 21 shows the pressure profile for different penalty factors. In n = 10E, the pressure profile is very smooth but the bumps which were discussed before appear in n = 100E and 1000E again. It is seen that the intensity of these bumps increase with the penalty factor. Another point is that we could solve the problem with t = 10−5 s for n = 10E and 100E but for n = 1000E we had to use t = 10−6 s. It is well known that increasing the penalty factor increases the condition number of the system and decreases the critical time step of the explicit methods. 4.3 A rectangular block with curved interface In order to show the capability of the method in the case of complicated interfaces, an elastic 2 m × 1 m block with a curved interface in the middle is considered; the interface’s curve is defined as y = 0.5 + 0.05 sin( 5π 2 x). The material properties are like the previous examples; the problem geometry and boundary conditions are shown in Fig. 22. The bottom face of the block is fixed clamped, the top face is fixed in y direction and it is subjected to a ramp displacement of 1m in 0.2 s in the x direction. The problem is in plane strain state, the normal contact factor is assumed to be

123

Comput Mech

n = 100E and the time step is taken to be t = 10−6 s. For the numerical simulations, four different mesh densities of 40 × 20, 60 × 30, 80 × 40 and 100 × 50 elements are considered. The problem geometry is complex and it involves large sliding and large strains. Figure 23 shows the σ yy stress contour for the four mesh densities at t = 0.1 s and t = 0.2 s. It is observed that the solutions are almost identical and they show convergence. High stress concentration occurs when the hills of the top and bottom interface are in contact. Figure 24, shows the normal pressure distribution on the upper face of the block at t = 0.1 s and t = 0.2 s. It is evident that as the interface moves, new regions come into contact and some regions become traction free. The results are in good agreement for different mesh densities. In order to investigate the performance of the formulation in very sever contact problems, the height of the interface is increased to 0.1 m thus the curve is defined as y = 0.5+0.1sin( 5π 2 x). The problem is solved for the 40×20 mesh and the boundary conditions are the same as before. Figure 25a shows the σ yy stress contour at t = 0.066 s; unlike the previous interface geometry, at some instances, the hanging node issue is observed in few elements; in fact, since

Fig. 20 Effect of penalty factor; two body contact with mesh 5 × 20 at t = 0.1 s a n = 10E, b n = 100E and c n = 1000E

123

the contact force is very high, even the node relocation strategy is not enough to suppress this phenomena. Figure 25b shows the same problem when node relocation is accompa-

Fig. 21 Effect of the penalty factor on pressure profile; two body contact with mesh 5 × 20 at t = 0.1 s

Comput Mech Fig. 22 The geometry and boundary conditions for a block with a curved interface problem

Fig. 23 the σ yy stress contour for a, c 40 × 20 mesh, c, d 60 × 30 mesh, e, f 80 × 40 mesh and g, h 100 × 50 mesh; left column t = 0.1 s and right column t = 0.2 s

nied with α artificial internal force (30); as it is observed the hanging node issue is fully removed but the stress field is changed a little in the corresponding elements. The contact pressure profile for both cases are compared in Fig. 26. The overall response is almost identical but the pressure in the elements which were prone to hanging node issue is changed slightly.

4.4 Impact of a ring to an obstacle The aim of the final example is to show the applicability of the proposed general framework to high speed impact problems. The problem is a 2D ring which is fired toward an elastic obstacle. The inner and outer radius of the ring are 0.3 and 0.4 m respectively and the obstacle has a dimen-

123

Comput Mech

Fig. 24 Contact pressure distribution for the four mesh densities a t = 0.1 s and b t = 0.2 s

Fig. 25 The σ yy stress contour for the block with a curve shape y = 0.5 + 0.1sin( 5π 2 x) interface at t = 0.066 s; a node relocation

is applied; hanging node problem is observed in this sever case and b node relocation is accompanied with α artificial internal force; the hanging node issue is removed but the stress field changes a little

sion of 1 m × 0.5 m. The ring has a Young’s modulus of E = 10 GPa, Poisson’s ratio of ν = 0.3 and density of ρ = 7800 kg/m3 . The corresponding values for the obstacle are E = 1000 GPa, ν = 0.3 and ρ = 7800 kg/m3 . The problem is in plane strain condition and Fig. 27 shows the geometry and mesh; It is seen that the meshes are structured and the contact area is defined by two interfaces, shown by gray lines; the first interface is placed close to the lower edge of the ring and the second one is close to the upper edge of the obstacle. The bottom face of the obstacle is fixed clamped and an initial vertical speed of v = −50 m/s is applied to the ring. The normal penalty factor is n = 100E and the time step is taken to be t = 10−7 s. Figure 28 shows the von-Mises stress contour at different time steps during the impact. When, the ring collides with

the obstacle, it deforms as the stress wave propagates throw the body and then it bounces back. It is obvious that the proposed formulation correctly captures the contact area and it applies suitable contact pressure to the interfaces. Note that the observable penetration of the ring in the obstacle, e.g. Fig. 28c is not the real penetration. As it was stated before, the contact is between the interfaces and not among the element edges. Thus, the penetrated parts of the elements are of no computational importance. In other words, in phantom node method, it is difficult to graphically show the results; since, we have the extra sections of the master and slave elements, which make the figures visually crowded. The same problem was solved for v0 = −100 m/s; Fig. 29 shows the corresponding von-Mises stress contour. It is seen that even in this high speed impact problem the contact for-

123

Comput Mech

Fig. 27 The geometry and mesh for the ring impact problem Fig. 26 The contact pressure profile for the block with a curve shape y = 0.5 + 0.1sin( 5π 2 x) interface at t = 0.066 s; comparison of the case in which only the node relocation remedy is used with the case for which node relocation is accompanied with α artificial internal force

mulation works properly. In impact problems conservation of energy plays a crucial role; Fig. 30 shows the time history of kinetic, internal and total energies of the system for

v0 = −50 m/s and v0 = −100 m/s. It is seen that as the ring collides with obstacle the kinetic energy is transformed into internal elastic energy and thereafter these energies are transformed into each other repeatedly; however, the total energy remains almost constant for both cases which shows

Fig. 28 von-Mises stress contour for initial speed of v0 = −50 m/s at a t = 0.003 s, b t = 0.006 s, c t = 0.009 s, d t = 0.012 s, e t = 0.0015 s and f t = 0.0018 s

123

Comput Mech

Fig. 29 von-Mises stress contour for initial speed of v0 = −100 m/s at a t = 0.003 s, b t = 0.006 s, c t = 0.009 s, d t = 0.012 s, e t = 0.0015 s and f t = 0.0018 s

Fig. 30 Energy time histories of kinetic, internal and total energies for a v0 = −50 m/s and b v0 = −100 m/s

123

Comput Mech

that energy is approximately conserved in this formulation. As it is expected the total energy in v0 = −100 m/s case is four times the v0 = −50 m/s case. More sophisticated formulations, like the mortar contact method proposed by Hesch and Betsch [30] can be combined with phantom-node approach to make sure of the perfect conservation of energy.

5 Conclusions In this paper a general frame work was presented for dynamic contact and impact problems based on the phantom-node extended finite element method. Large sliding was implemented with a master-slave approach and large deformations and large strains were considered based on an updated Lagrangian scheme. The temporal discretization was performed with central difference method. Implementational details and obstacles were discussed. The proposed framework removes the need for complicated meshing and the interfaces are separated from the geometry. The method is simple to implement in existing finite element codes and the extension to 3D problems is almost straightforward. Several numerical examples were presented which showed the application of the method in different complex contact-impact problems.

References 1. Kobayashi S, Kobayashi S, Oh S-I, Altan T (1989) Metal forming and the finite-element method, vol 4. Oxford University Press on Demand, Oxford 2. Ambrosio JAC (2003) Contact and impact models for vehicle crashworthiness simulation. Int J Crashworthiness 8(1):73–86 3. Lim CT, Shim VPW, Ng YH (2003) Finite-element modeling of the ballistic impact of fabric armor. Int J Impact Eng 28(1):13–31 4. Yao H, Gao H (2006) Mechanics of robust and releasable adhesion in biology: bottom-up designed hierarchical structures of gecko. J Mech Phys Solids 54(6):1120–1146 5. Liu F, Borja RI (2009) An extended finite element framework for slow-rate frictional faulting with bulk plasticity and variable friction. Int J Numer Anal Meth Geomech 33(13):1535–1560 6. Dolbow J, Moës N, Belytschko T (2001) An extended finite element method for modeling crack growth with frictional contact. Comput Methods Appl Mech Eng 190(51):6825–6846 7. Hughes TJR, Cottrell JA, Bazilevs Y (2005) Isogeometric analysis: cad, finite elements, nurbs, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 194(39):4135–4195 8. Belytschko T, Black T (1999) Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 45(5):601–620 9. Khoei AR, Shamloo A, Azami AR (2006) Extended finite element method in plasticity forming of powder compaction with contact friction. Int J Solids Struct 43(18):5421–5448 10. Khoei AR, Nikbakht M (2007) An enriched finite element algorithm for numerical computation of contact friction problems. Int J Mech Sci 49(2):183–199

11. Elguedj T, Gravouil A, Combescure A (2007) A mixed augmented lagrangian-extended finite element method for modelling elasticplastic fatigue crack growth with unilateral contact. Int J Numer Meth Eng 71(13):1569–1597 12. Liu F, Borja RI (2008) A contact algorithm for frictional crack propagation with the extended finite element method. Int J Numer Methods Eng 76(10):1489–1512 13. Liu F, Borja RI (2010) Stabilized low-order finite elements for frictional contact with the extended finite element method. Comput Methods Appl Mech Eng 199(37):2456–2471 14. Hirmand M, Vahab M, Khoei AR (2015) An augmented lagrangian contact formulation for frictional discontinuities with the extended finite element method. Finite Elem Anal Des 107:28–43 15. Geniaut S, Massin P, Moës N (2007) A stable 3D contact formulation using X-FEM. Eur J Comput Mech 16(2):259–275 16. Nistor I, Guiton MLE, Massin P, Moës N, Geniaut S (2009) An X-FEM approach for large sliding contact along discontinuities. Int J Numer Methods Eng 78(12):1407–1435 17. Liu F, Borja RI (2010) Finite deformation formulation for embedded frictional crack with the extended finite element method. Int J Numer Methods Eng 82(6):773 18. Khoei AR, Mousavi SMT (2010) Modeling of large deformationlarge sliding contact via the penalty X-FEM technique. Comput Mater Sci 48(3):471–480 19. Song J-H, Areias P, Belytschko T (2006) A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer Methods Eng 67(6):868–893 20. Nistor I, Pantalé O, Caperaa S (2008) Numerical implementation of the extended finite element method for dynamic crack analysis. Adv Eng Softw 39(7):573–587 21. Chau-Dinh T, Zi G, Lee P-S, Rabczuk T, Song J-H (2012) Phantomnode method for shell models with arbitrary cracks. Comput Struct 92:242–256 22. de Souza Neto EA, Peric D, Owen DRJ (2011) Computational methods for plasticity: theory and applications. Wiley, Hoboken 23. Khoei AR (2014) Extended finite element method: theory and applications. wiley 24. Menouillard T, Rethore J, Moes N, Combescure A, Bung H (2008) Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation. Int J Numer Methods Eng 74(3):447–474 25. de Borst R, Remmers JJC, Needleman A (2006) Mesh-independent discrete numerical representations of cohesive-zone models. Eng Fract Mech 73(2):160–177 26. Hansbo A, Hansbo P (2004) A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Methods Appl Mech Eng 193(33–35):3523–3540 27. Rabczuk T, Belytschko T (2007) A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 196(29):2777–2799 28. Talebi H, Samaniego C, Samaniego E, Rabczuk T (2012) On the numerical stability and mass-lumping schemes for explicit enriched meshfree methods. Int J Numer Methods Eng 89(8):1009– 1027 29. Flanagan DP, Belytschko T (1981) A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. Int J Numer Methods Eng 17(5):679–706 30. Hesch C, Betsch P (2011) Transient three-dimensional contact problems: mortar method. mixed methods and conserving integration. Comput Mech 48(4):461–475

123