INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2000; 00:1–6 Prepared using nmeauth.cls [Version: 2002/09/18 v2.02]
A meshfree thin shell method for nonlinear dynamic fracture
T.Rabczuk
†
, P.M.A.Areias
+
, T.Belytschko
∗,k
Department of Mechanical Engineering, Northwestern University, Evanston, IL 60208-311, U.S.A
SUMMARY
A meshfree method for thin shells with finite strains and arbitrary evolving cracks is described. The
C
1
displacement continuity requirement is met by the approximation, so no special treatments for
fulfilling the Kirchhoff condition are necessary. Membrane locking is eliminated by the use of a cubic
or quartic polynomial basis. The shell is tested for several elastic and elasto-plastic examples and
shows good results. The shell is subsequently extended to modelling cracks. Since no discretization
of the director field is needed, the incorporation of discontinuities is easy to implement and straight
forward. Copyright
c
° 2000 John Wiley & Sons, Ltd.
key words: meshfree methods, cracks, cohesive models, KL constraint, shell
1. INTRODUCTION
This paper describes a meshfree method for treating the dynamic fracture of shells. It includes
both geometric and material nonlinearities and also includes a meshfree fluid model, so that
complex fluid-structure interaction problems are feasible. Here we illustrate this capability
with the fracture of a fluid-filled cylinder that is impacted by a penetrating projectile. The
shell formulation is based on the Kirchhoff-Love (KL) theory.
A meshfree thin shell based on the imposition of the KL constraints was first proposed by
Krysl and Belytschko [1]. However, the shell was developed for small deformations, small strains
and elasticity. Three dimensional modelling of shear deformable shells and degenerated shells
in a meshfree context was studied by Noguchi et al. [2], Li et al. [3] and Kim et al. [4]. Usually,
a low order polynomial basis was used, e.g. in [3], a trilinear polynomial basis was applied
and the method was applied to several non-linear problems. However, for thin structures,
three dimensional modelling of shell structures is computationally expensive. Garcia et al. [5]
developed meshfree methods for plates and beams; the higher continuity of meshfree shape
functions was exploited for Mindlin-Reisner plates in combination with a p-enrichment. Wang
and Chen [6] proposed a meshfree method for Mindlin-Reisner plates. Locking is treated using
∗
Correspondence to: T. Belytschko, Department of Mechanical Engineering, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208-311, U.S.A. E-mail: t-belytschko@nwu.edu
k Walter P.Murphy, Professor of Computational Mechanics
†, + Post-Doctoral Research Fellow, Department of Mechanical Engineering
Received
Copyright
c
° 2000 John Wiley & Sons, Ltd. Revised
2 T. RABCZUK, P.M.A. AREIAS, T. BELYTSCHKO
second order polynomials for the approximation of the translational and rotational motion in
combination with a curvature smoothing stabilization. Kanok-Nukulchai et al. [7] addressed
shear locking for plates and beams in the element-free Galerkin method.
We develop a meshfree thin shell that combines classical shell theory with a continuum based
shell. The kinematic assumptions of classical KL shell theory is adopted. We make use of the
generality provided by the continuum description, so that constitutive models developed for
continua are easily applicable to shells. The formulation is valid for finite strains.
We include in the shell the capability to model cracks, which are modeled either by cracked
particles as in Rabczuk and Belytschko [8] or by a local partition of unity [9, 10, 11]. Due
to the higher order continuity of meshfree methods, which enables the use of Kirchhoff-
Love shell theories in pristine form, the incorporation of discontinuities is very simple and
straight forward. The director field is not discretized, which simplifies the incorporation of
discontinuities. In our meshfree model, the concepts for modelling cracks in continua can be
adopted directly to shells.
The paper is arranged as follows: First, the kinematics of the shell is described. Then, the
meshfree method, the element-free Galerkin (EFG) method, is reviewed and the concept how
to incorporate continuum constitutive models is described. The extensions to modelling cracks
is described in section 5. Finally, we test the meshfree shell for different elastic linear and
nonlinear problems, plastic problems and cracking problems.
2. SHELL MODEL
2.1. Kinematics
Consider a body Ω with material points X ∈ Ω
0
of the shell in the reference configuration
with discontinuities, e.g. cracks, on lines Γ
c
0
. The boundary is denoted by Γ
0
, where Γ
u
0
and
Γ
t
0
are the complementary boundaries on which, respectively, displacements and tractions are
prescribed. We consider a surface parametrized by two independent variables θ
α
, α = 1 to 2;
the surface in the reference (initial) configuration is described by R(θ
α
), R ∈ <
3
. The material
points in the reference configuration are given by
X(θ
i
) = R(θ
α
) + θ
3
d
2
N(θ
α
) (1)
where θ
i
are curvilinear coordinates, −1 ≤ θ
3
≤ 1, d is the thickness of the shell, N is the
shell normal and R is a point on the mid-surface of the shell in the reference configuration S
0
.
Upper Latin and Greek indices range from 1 to 3 and from 1 to 2, respectively and refer to
quantities in the cartesian or curvilinear coordinate system. The current configuration is given
by
x(θ
i
, t) = r(θ
α
, t) + θ
3
d
2
n (θ
α
, t) (2)
where t is the time, n is the director field and r is a point on the mid-surface position in the
current configuration. The first and second fundamental forms are given by
A
αβ
= R
,α
· R
,β
(3)
B
αβ
= R
,α,β
· N = −R
,α
· N
,β
(4)
Copyright
c
° 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1–6
Prepared using nmeauth.cls
A DEMONSTRATION OF THE INT. J. NUMER. METH. ENGNG CLASS FILE WITH FLUID-STRUCTURE INTERACTION3
The curvilinear coordinates θ
α
are such that R
,α
form a basis for the tangent space in X ∈ S
0
.
For arbitrary θ
3
, we define a family of surfaces S(θ
3
) with S
0
= S(0) for which the tangent basis
is established from (1) as X
,α
= R
,α
+
θ
3
2
N
α
. We extend this basis by including N = X
,3
2
d
.
The resulting basis spans Ω
0
, and we can then define the metric of Ω
0
as G
ij
= X
,i
· X
,j
. The
dual basis is given by G
i
= G
ij
X
,j
with
£
G
ij
¤
= [G
ij
]
−1
.
The Cauchy-Green tensor is
C = F
T
F = (x
,i
· x
,j
)
| {z }
C
ij
G
i
⊗ G
j
(5)
where x are the spatial position coordinates and F is the deformation gradient.
The Kirchhoff-Love hypothesis is imposed by requiring that n is perpendicular to r
,α
,
α = 1, 2:
n =
r
,1
× r
,2
kr
,1
× r
,2
k
(6)
2.2. Virtual work
The weak form of the momentum equation is written with the principle of virtual work (see
e.g. Belytschko et al. [12]): find r ∈ V such that
δW = δW
int
− δW
ext
+ δW
kin
− δW
E
= 0 ∀δr ∈ V
0
(7)
where
V =
©
r(·, t)|r(·, t) ∈ H
2
(Ω
0
/Γ
c
0
), r(·, t) =
¯
r(t) on Γ
u
0
, r discontinuous on Γ
c
0
ª
V
0
= {δr|δr ∈ V, δr = 0 on Γ
u
0
, δr discontinuous on Γ
c
0
} (8)
δW
int
=
Z
Ω
0
\Γ
c
0
©
s
αβ
x
,α
· δx
,β
[G
1
· (G
2
× G
3
)]
ª
dΩ (9)
δW
ext
=
Z
Ω
0
\Γ
c
0
%
0
b · δu dΩ
0
+
Z
Γ
t
0
¯
t
0
· δu dΓ
0
(10)
δW
E
=
Z
Γ
c
¯
t
c
· δ[[u]] dΓ (11)
δW
kin
=
Z
Ω
0
\Γ
c
0
%
0
δu ·
¨
u dΩ
0
(12)
where the prefix δ identifies the test function and W
ext
is the external energy, W
int
designates
the internal energy, W
E
is the crack cohesive energy and W
kin
the kinetic energy, %
0
is the
density, s is the Kirchhoff stress, b is the body force and
¯
t
0
the prescribed traction; superposed
dots denote material time derivatives.
For details of the discrete equations see Rabczuk and Belytschko [8, 13, 14], Belytschko et al.
[15].
Copyright
c
° 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1–6
Prepared using nmeauth.cls
4 T. RABCZUK, P.M.A. AREIAS, T. BELYTSCHKO
2.3. Discretization
The approximation of the shell surface is given by
r(θ
α
, t) =
X
I∈W
Φ
I
(θ
α
) r
I
(t) (13)
where Φ
I
(θ
α
) are the shape functions and W is the domain of influence of the corresponding
particle.
We require C
3
displacement continuity in the meshfree method. So there is no need to
discretize the director field n and n is readily obtained from eq. (6). The variation of the
motion x (and their spatial derivatives) is given by
δx = δr + θ
3
d
2
δn (14)
The variation of the normal can be expressed in terms of r:
δn = a
1
a
4
· a
2
(15)
where a
1
=
1
kr
,1
×r
,2
k
, a
2
= r
,1
× δr
,2
− r
,2
× δr
,1
and a
4
= (I − n ⊗ n).
The derivatives of the motion x are given by
x
,α
=
X
I∈W
Φ
I,α
r
I
|
{z }
r
,α
+θ
3
d
2
n
,α
(16)
with
n
,γ
= a
1
a
4
· a
5
(17)
with a
5
= r
,1
× r
,2γ
− r
,2
× r
,1γ
.
To obtain the variation δx
,α
, the derivatives of the variation of the normal have to be
computed. Considering eq. (15), the derivatives are:
δn
,γ
= − a
3
1
[(r
,1
× r
,2
) · a
5
](a
4
· a
2
) − 2a
1
(n ⊗ a
5
)
S
· a
2
+
a
1
a
4
· δa
5
(18)
3. EFG-MESHFREE METHOD
We use the element free Galerkin (EFG) method that is based on a moving least square (MLS)
approximation, see e.g. Belytschko et al. [16], Belytschko and Lu [17]. For the meshfree shell, at
least a quadratic polynomial basis is required. Krysl and Belytschko [1] showed that a quadratic
basis can lead to membrane locking in meshfree thin shells and showed that a quartic basis
completely eliminates membrane locking. A quartic polynomial basis can be written as:
p(θ
α
) =
¡
1, θ
1
, θ
2
, (θ
1
)
2
, θ
1
θ
2
, (θ
2
)
2
, (θ
1
)
3
, (θ
1
)
2
θ
2
, (θ
2
)
2
θ
1
, ...
... (θ
2
)
3
, (θ
1
)
4
, (θ
1
)
3
θ
2
, (θ
1
)
2
(θ
2
)
2
, θ
1
(θ
2
)
3
, (θ
2
)
4
¢
(19)
Copyright
c
° 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1–6
Prepared using nmeauth.cls
A DEMONSTRATION OF THE INT. J. NUMER. METH. ENGNG CLASS FILE WITH FLUID-STRUCTURE INTERACTION5
For our examples, it was sufficient to use a cubic basis.
The MLS approximation is given by:
r
h
i
(θ
α
, t) =
X
J∈S
p(θ
α
J
) a
i
(θ
α
, t) (20)
with a
i
(θ
α
, t) chosen to minimize
J
i
=
X
J∈S
¡
p(θ
α
J
)
T
a
i
(θ
α
, t) − r
Ji
(t)
¢
2
w(θ
α
− θ
α
J
, h). (21)
with respect to a
i
(θ
α
, t). This leads to the approximation
r
h
i
(θ
α
, t) =
X
J∈S
r
Ji
(t) Φ
J
(θ
α
) (22)
with
Φ
J
= p
T
(θ
α
) · M
−1
(θ
α
) · p(θ
α
J
) w(θ
α
− θ
α
J
, h) (23)
M(θ
α
) =
X
J∈S
p(θ
α
J
) p
T
(θ
α
J
) w(θ
α
− θ
α
J
, h). (24)
w(r) = w(θ
α
J
−θ
α
I
, h) is the kernel function that determines the order of continuity. In the EFG
method, the continuity of the shape functions is equivalent to the continuity of the kernel, see
[18]. A cubic spline kernel leads to C
2
continuity, a quartic spline to C
3
continuity. Recall that
for a Kirchhoff shell, C
1
continuity is required.
For rectangular plates and for cylinders, the Jacobian is constant and hence θ
i
are linear
combinations of x and linear completeness is guaranteed. It is noted that the same shape
functions are employed in shape and displacement approximations to guarantee strain-free
states in rigid body motion, see Krysl and Belytschko [1].
4. CONTINUUM CONSTITUTIVE MODELS
For the constitutive model, we adopt the algorithm of table I. We use a two-dimensional non-
symmetric radial return and rotate so that the 3-3 component corresponds to the normal. The
normal strain (and consequently the normal stress) is filtered, according to figure 1. For more
details on that specific constitutive model, the reader is referred to Areias et al. [19].
5. CRACK MODELING
The approach is extended to shells with cracks by enriching the mid-surface motion r(θ
α
, t)
with a discontinuous function. The jump in the director field is then obtained directly via the
discontinuous part of r(θ
α
, t). The fact that there is no need to discretize the director field,
facilitates the incorporation of discontinuities in shells.
The force introduced here corresponds to the resistance to opening, which is a function
of the opening displacement. The opening displacement can be written as a function of the
mid-surface position and the director.
Copyright
c
° 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1–6
Prepared using nmeauth.cls
![Figure 8. Load-displacement curve of the pinched clamped cylinder compared with Areias et al. [26], Ibrahimbegović et al. [34], Eriksson and Pacoste [35]](/figures/figure-8-load-displacement-curve-of-the-pinched-clamped-2az09ya9.webp)
![Figure 9. Load-displacement curve of the open pullout cylinder compared with Areias et al. [26], Peng and Crisfield [38], Masud et al. [39]](/figures/figure-9-load-displacement-curve-of-the-open-pullout-2mmha705.webp)
