A Meshing Strategy for a Quadratic Iso-parametric FEM in Cavitation Computation in Nonlinear Elasticity

0.44 MB·

by  Chunmei Su

A Meshing Strategy for a Quadratic
7
1 Iso-parametric FEM in Cavitation Computation
0
2
n in Nonlinear Elasticity∗
a
J

5

]
A Chunmei Su and Zhiping Li†
N
.
h
t
a
m
LMAM & School of Mathematical Sciences,
[

1
v Peking University, Beijing 100871, China
3
3
2
1
0
. Abstract
1
0
The approximation properties of a quadratic iso-parametric finite ele-
7
1
ment method for a typical cavitation problem in nonlinear elasticity are
:
v
analyzed. More precisely, (1) the finite element interpolation errors are es-
i
X
tablishedintermsofthemeshparameters; (2)ameshdistributionstrategy
r
a
based on an error equi-distribution principle is given; (3) the convergence
of finite element cavity solutions is proved. Numerical experiments show
that,infact,theoptimalconvergenceratecanbeachievedbythenumerical
cavity solutions.
Keywords: cavitation, quadratic iso-parametric FEM, error analysis, meshing strat-
egy, convergence, nonlinear elasticity.
AMS Subject Classification: 65N12, 65N15, 65N30, 65N50, 74B20, 74G15,74M99.
∗The researchwas supported by the NSFC project 11171008and RFDP of China.
†Corresponding author, email: [email protected]
1
1 Introduction
Nonlinear soft elastic materials, such as polymers, biological tissues, rubbers,
etc., can display a particular singular deformation, referred to in the literature as
cavitation, when strong external force is applied [6, 7, 9, 17, 25]. The occurrence
and growth of cavities is considered closely related to the material instability and
to the damage and failure mechanisms of the materials [10, 14, 15, 16, 19]. A
huge number of work has been done by numerous authors analysing cavitation
experimentally, theoretically as well as numerically.
Generally speaking, there are two representative approaches characterizing
cavitation. One is the so-called defect model, which is based on the hypothesis
that cavities grow from pre-existing micro defects. Under this assumption, Gent
and Lindley [9] analysed the critical hydrostatic pressure at which a given unit
spherical void in an infinite extension of a Neo-Hookean material would blow up,
which was in a good agreement with their experiments therein. The other is the
perfect model established by Ball [3] based on the analytical evidence that, under
certain circumstances, a deformation with cavities created in an originally intact
material can be energetically favorable. It is shown that, under the assumption
that the cavities can appear only at a finite number of fixed points in the intact
materials, the solution of the defect model converges to the solution of the perfect
model [11, 32] as the radii of the pre-existing small voids go to zero. In addition,
analytical and numerical evidences indicate that whether a point can serve as a
possible position of cavitation can be evaluated by calculating the corresponding
configurational forces [22, 31].
The perfect model typically exhibits the Lavrentiev phenomenon [18] when
there is a cavitation solution, leading to the failure of the conventional finite
element methods [1, 4]. Though there are existing numerical methods developed
to deal with the Lavrentiev phenomenon [1, 4, 20, 28], they do not seem to be
suitable to tackle the cavitation problem. In fact, most of the numerical studies
on cavitation are based on the defect model, in which one considers to minimize
2
the total energy of the form
E(u) = W(∇u(x))dx, (1.1)
Z
Ω̺
in the set of admissible functions
A = {u ∈ W1,p(Ω ;Rn) is one-to-one a.e. : u| = u ,det∇u > 0 a.e.}, (1.2)
̺ Γ0 0
where Ω = Ω\ m B (a ) ⊂ Rn(n = 2,3) denotes the region occupied by an
̺ i=1 ̺i i
elastic body in iSts reference configuration, B (a ) = {x ∈ Rn : |x − a | < ̺ }
̺i i i i
is the pre-existing spherical hole centered at a with small radius ̺ > 0, W :
i i
Mn n → R+ is the elastic stored energy density of the material, Mn n denotes
+× +×
the set of n×n matrices with positive determinant, Γ is the boundary of Ω.
0
A typical example of the elastic stored energy density is of the form
W(F) = ω|F|p+g(detF), ∀F ∈ Mn n, (1.3)
+×
where ω > 0 is a material constant, n−1 < p < n, and g : (0,∞) → (0,∞) is a
continuously differentiable strictly convex function characterizing the compress-
ibility of the material and satisfies
g(d)
g(d) → +∞ as d → 0, and → +∞ as d → +∞. (1.4)
d
As was shown by Ball [3], this kind of functional can have a singular minimizer
displaying cavitation. Further studies on the existence of singular minimizers in
Sobolev spaces are referred to [11, 27, 30].
One of the main difficulties in the computation of immense growth of voids
is the orientation-preservation of the finite element deformation, which is a cru-
cial constraint and is characterized by the pointwise positivity of the Jacobian
determinant of the deformation gradient. For the conforming piecewise affine
finite element, the condition leads to an unbearably large amount of degrees
of freedom [35]. Some numerical methods [21, 22, 23, 35] have been designed
to overcome the difficulty, and have shown some success numerically. However,
strict analytical results are insufficient. The only practical analytical results for
3
the cavitation computation known to the authors so far are [34], where a suf-
ficient orientation-preservation condition and the interpolation error estimates
were given for a dual-parametric bi-quadratic finite element method, and [33],
where a set of sufficient and necessary orientation-preservation conditions for the
quadratic iso-parametric finite element interpolation functions of radially sym-
metric cavity deformations are derived.
In this paper, we analyze the approximation properties of a quadratic iso-
parametric finite element for the typical cavitation problem. The analytical re-
sults on the errors of finite element interpolation functions lead to a delicate
relationship between the elastic energy error and the mesh parameters, which
together with the orientation-preservation conditions (see Remark 3.4 and [33])
enable us to establish a mesh distribution strategy guaranteeing that the cor-
responding finite element cavitation solution is orientation preserving and its
relative error on the elastic energy is O(h2), where h is the mesh size in the far
field, i.e. a given distance away from the cavity. Above all, for the first time to
our knowledge, the convergence of the finite element cavitation solutions in W1,p
norm is proved. In fact, the numerical experiments show that the optimal order
of convergence rate is achieved by the numerical cavitation solutions obtained on
the meshes produced by our meshing strategy.
Since the cavitation solution is generally considered to be quite regular except
in a neighborhood of the voids, where the material experiences extremely large
expansion dominant deformations and the difficulty of the computation as well as
analysis lies, we restrict ourselves to a simplified problem with Ω = B (0)\B (0)
̺ 1 ̺
in R2 and a simple expansionary boundary condition given by u = λx.
0
Thestructureofthepaperisasfollows. In§2, weintroducetheiso-parametric
finite element method and a radially symmetric large expansion accommodating
triangulation method briefly. § 3 is devoted to analyzing the interpolation errors
of the cavitation solutions. The meshing strategy is given in § 4, where the
convergence theorem is also established. The numerical results are presented in
§ 5. Some concluding remarks are made in § 6.
4
2 Preliminaries
2.1 The quadratic iso-parametric FEM
Let (Tˆ,Pˆ,Σˆ) be a quadratic Lagrange reference element. Define F : Tˆ → R2
T
F ∈ Pˆ2 = (P (Tˆ))2,
T 2

 3 (2.1)

x = F (xˆ) = a µˆ (xˆ)+ a µˆ (xˆ),
T i i ij ij
X X
i=1 1 i<j 3
 ≤ ≤

where a ,1 ≤ i ≤ 3, and a ,1 ≤ i < j ≤ 3 are given points in R2, and
i ij
µˆ (xˆ) = λˆ (xˆ)(2λˆ (xˆ)−1), µˆ (xˆ) = 4λˆ (xˆ)λˆ (xˆ),
i i i ij i j
with λˆ (xˆ),1 ≤ i ≤ 3 being the barycentric coordinates of Tˆ. If the map F
i T
defined above is an injection, then, T = F (Tˆ) is a curved triangular element
T
as shown in Figure 2. The standard quadratic iso-parametric finite element is
defined as a finite element triple (T,P ,Σ ) with
T T
T = F (Tˆ) being a curved triangle element,
T

PT = {p : T → R | p = pˆ◦FT−1, pˆ∈ Pˆ = P2(Tˆ)}, (2.2)
Σ = {p(a ),1 ≤ i ≤ 3;p(a ),1 ≤ i < j ≤ 3}.
T i ij




xˆ2 x2
a a
aˆ3 a 13 3
1 1
a a
12 23
aˆ23
aˆ13
a
aˆ1 aˆ12 aˆ2 2
Oˆ 1 xˆ1 O x1
Figure 1: The reference element Tˆ. Figure2: Acurved triangularelement T.
5
2.2 Large expansion accommodating triangulations
Let’s have a look at how the iso-parametric FEM is applied in cavitation com-
putation. As introduced in [22], let x (k = 1,··· ,m) be the center of the k−th
k
defect with small radius ̺ . In the domain far away from the defects where the
k
deformation is regular, the mesh is given by general straight edged triangulation.
To accommodate the large expansionary deformation around the defects, the tri-
angulation near the defects is given by the quadratic mapping as in (2.1), where
given 3 vertices a , denote (r (x),θ (x)) the local polar coordinates of x with
i k k
respect to x , set
k
a = (r cosθ ,r sinθ )+x , (2.3)
ij ij ij ij ij k
where
r(a )+r(a ) θ(a )+θ(a )
i j i j
r = , θ = .
ij ij
2 2
Moreover, the elements on the boundary are also adjusted to achieve a better
approximation of the region. With this kind of curved elements near the defects
and general straight triangles elsewhere, the mesh can better accommodate the
locally large expansionary deformations. As an example, an EasyMesh produced
mesh J with m = 2 is shown in Figure 3, and the final mesh J (see Figure 4)
′
is obtained by adding to J around each defect two layers of radially symmetric
′
mesh of the kind shown in Figure 5(a), which is a standard curved triangulation
around a prescribed circular ring with inner radius ǫ = 0.01 and thickness τ =
0.01. Figure 5(b) shows that an outer layer of standard curved triangulation is
connected to a doubly refined inner one by a conforming layer of nonstandard
curved triangulation. For the convenience of reference, we classify the curved
triangular elements in Figure 5 into four basic types, and denote them as types
A, B, C and D.
3 Interpolation errors of cavity deformations
There are standard error estimates [5] on the interpolation functions in the
iso-parametric finite element function spaces, however, because of the highly
6
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1 −1
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
Figure 3: An EasyMesh J . Figure 4: J, a mesh adapted to cavity.
′
0.02 0.04
0.015 0.03
0.01 0.02
0.005 0.01 C
A
B
0 ε τ 0 AB D
−0.005 −0.01
−0.01 −0.02
−0.015 −0.03
−0.02 −0.04
−0.015 −0.01 −0.005 0 0.005 0.01 0.015 −0.04 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04
(a) A standard layer of curved triangulation. (b) A conforming layer links 2 standard ones.
Figure 5: Triangulation layers with elements of types A, B, C and D.
anisotropic deformation inevitably involved in the cavitation computation, the
error bounds so obtained depend generally on the size of the initial void ǫ and
0
blow up when ǫ → 0. In this section, with some more sophisticated manipula-
0
tions and calculations, we are able to establish the ǫ -independent interpolation
0
error estimates, including those onthedeformation anditsJacobiandeterminant,
and in particular the elastic energy, for the radially symmetric cavity deforma-
7
tions in the quadratic iso-parametric finite element function spaces defined on the
meshes consisting of only elements of types A, B, C and D.
To have a better picture in our mind for the problem given below, we first
introduce some notations. Let ǫ and τ represent respectively the inner radius
and the thickness of the circular annulus as shown in Figure 5(a), let N be the
number of the evenly spaced nodes on the outer circle of the circular ring, and
denote Ω = {x ∈ R2 : ǫ ≤ |x| ≤ ǫ+τ}. Throughout the paper, the notation
(ǫ,τ)
Φ . Ψ means that there exists a generic constant C independent of ǫ and τ such
that |Φ| ≤ CΨ, and Φ ∼ Ψ means that Ψ . Φ . Ψ.
3.1 The interpolation function and its Jacobian
x
For a radially symmetric function v(x) = s(|x|)x, the iso-parametric finite ele-
| |
ment interpolation function can be written as (see § 2)
3
Πv(x) = b µˆ (xˆ)+ b µˆ (xˆ), (3.1)
i i ij ij
X X
i=1 1 i<j 3
≤ ≤
wherexˆ = F 1(x),andwhere, forarepresentative oftypeAelement, b = (s ,0),
T− 1 0
b = s (cos π,−sin π), b = s (cos π,sin π), b = s (cos π ,−sin π ), b =
2 1 N N 3 1 N N 12 1/2 2N 2N 13
s (cos π ,sin π ), b = (s ,0). On this element, let y = xˆ +xˆ and z = xˆ xˆ ,
1/2 2N 2N 23 1 1 2 1 2
we have
π
Πv(x) = s +α y +2α y2 −4s sin2 (xˆ2 +xˆ2),(2γy −β)(xˆ −xˆ ) , (3.2)
(cid:16) 0 2 1 1 2N 1 2 2 1 (cid:17)
where
π
α = s +s −2s cos , (3.3)
1 0 1 1/2
2N
π π
α = −3s −s cos +4s cos , (3.4)
2 0 1 1/2
N 2N
π π
β = s sin −4s sin , (3.5)
1 1/2
N 2N
π π
γ = s sin −2s sin . (3.6)
1 1/2
N 2N
8
Hence
∂Πv  α−8s1xˆ1sin2 2πN α−8s1xˆ2sin2 2πN 
= , (3.7)
∂xˆ
 β −4γxˆ −β +4γxˆ 
1 2
 
where α = 4α y +α . It follows that
1 2
∂Πv π
det (xˆ ,xˆ ) = H(y,z) , 16γα y2 −64s γsin2 z
∂xˆ 1 2 1 1 2N
π
+ −8β α −s sin2 +4γα y −2βα . (3.8)
1 1 2 2
2N
(cid:0) (cid:0) (cid:1) (cid:1)
On a representative element of type B with b = (s ,0), b = s (cos π,sin π),
1 1 2 0 N N
b = s (cos π,−sin π), b = s (cos π ,sin π ), b = s (cos π ,−sin π ),
3 0 N N 12 1/2 2N 2N 13 1/2 2N 2N
b = (s ,0), again denote y = xˆ +xˆ and z = xˆ xˆ , one has
23 0 1 2 1 2
π
Πv(x) = s +α¯ y +2α¯ y2 −4s sin2 (xˆ2 +xˆ2),(2γ¯y −β¯)(xˆ −xˆ ) , (3.9)
(cid:16) 1 2 1 0 2N 1 2 1 2 (cid:17)
where
π
α¯ = s +s −2s cos , (3.10)
1 0 1 1/2
2N
π π
α¯ = −3s −s cos +4s cos , (3.11)
2 1 0 1/2
N 2N
π π
¯
β = s sin −4s sin , (3.12)
0 1/2
N 2N
π π
γ¯ = s sin −2s sin . (3.13)
0 1/2
N 2N
Hence
∂Πv  α¯ −8s0xˆ1sin2 2πN α¯ −8s0xˆ2sin2 2πN 
= , (3.14)
∂xˆ
¯ ¯
 −β +4γ¯xˆ β −4γ¯xˆ 
1 2
 
∂Πv π
det = H(y,z) = −16γ¯α¯ y2+64s γ¯sin2 z
∂xˆ 1 0 2N
π
+(8β¯(α¯ −s sin2 )−4γ¯α¯ )y +2β¯α¯ . (3.15)
1 0 2 2
2N
On a representative element of type C with b = (s ,0), b = (s ,0), b =
1 0 2 1 3
s (cos π,sin π),b = (s ,0),b = s (cos π ,sin π ),b = s (cos π ,sin π ),
0 N N 12 1/2 13 0 2N 2N 23 1/2 2N 2N
9
one has
π π
Πv(x) = s +α˜ xˆ +s α˜ xˆ +2α˜ xˆ2−8sin2 xˆ s cos xˆ −(s −s )xˆ ,
(cid:16) 0 1 1 0 2 2 3 1 4N 2 0 2N 2 0 1/2 1
(cid:0) (cid:1)
π π π
2sin xˆ s (2−cos )−4s sin2 xˆ +2(s −s )xˆ , (3.16)
2 0 0 2 1/2 0 1
2N 2N 4N (cid:17)
(cid:0) (cid:1)
where α˜ = 4s −s −3s , α˜ = 4cos π −cos π −3, α˜ = s +s −2s .
1 1/2 1 0 2 2N N 3 0 1 1/2
Similarly, on a representative element of type D with b = (s ,0), b =
1 0 2
s (cos π,−sin π), b = (s ,0), b = s (cos π ,−sin π ), b = (s ,0), b =
0 N N 3 1 12 0 2N 2N 13 1/2 23
s (cos π ,−sin π ), one has
1/2 2N 2N
π π
Πv(x) = s +s α˜ xˆ +α˜ xˆ +2α˜ xˆ2−8sin2 xˆ s cos xˆ −(s −s )xˆ ,
(cid:16) 0 0 2 1 1 2 3 2 4N 1 0 2N 1 0 1/2 2
(cid:0) (cid:1)
π π π
−2sin xˆ s (2−cos )−4s sin2 xˆ +2(s −s )xˆ . (3.17)
1 0 0 1 1/2 0 2
2N 2N 4N (cid:17)
(cid:0) (cid:1)
x
Throughout this section, we assume that u(x) = r(|x|)x represent a cavity
| |
deformation, where r(R) defined on (0,1] is smooth, positive, increasing, and
convex with r(j)(R)(j = 0,1,2,3) bounded. Moreover, it satisfies inf r(R) > 0,
R (0,1]
∈
mR ≤ r (R) ≤ MR, with 0 < m < M which is shown to be naturally satisfied for
′
the energy minimizer of (1.3) in [33]. For simplicity of the notations, we denote
y = xˆ +xˆ ∈ [0,1], z = xˆ xˆ ∈ [0, y2], where xˆ = (xˆ ,xˆ ) ∈ Tˆ.
1 2 1 2 4 1 2
3.2 The error of the interpolation function
x
We will estimate, in this subsection, the errors between u(x) = r(|x|)x and
| |
its interpolation function Πu(x) in the polar coordinates. Let (R,θ) be the
polar coordinates of x on the reference configuration and let (ζ,ϕ) be the polar
x
coordinates of the interpolation function Πu(x) of u(x) = r(|x|)x = (r(R),θ) on
| |
the deformed configuration.
Lemma 3.1 Denote x = F (xˆ) = R(cosθ,sinθ), then for the typical element A,
T
one has
R = ǫ+τy +O(τN 2 +ǫN 4) = (ǫ+τy)(1+O(N 2)), (3.18)
− − −
π
θ = (xˆ −xˆ )+O(N 3), (3.19)
2 1 −
N
π π ǫ
det∇x = 4τ sin (ǫ+τy +2τ sin2 ) 1+O(N 2 + N 4) ; (3.20)
− −
2N 4N τ
(cid:0) (cid:1)
10