Stabilized Finite Element Formulations for Shells in a Bending Dominated State.
We consider the design of finite element methods for the Naghdi shell model in the case when the deformation is  bending dominated. Two formulations based on  stabilizing techniques are introduced and it is proved that they are stable, hence free from locking. The theoretical estimates are confirmed by numerical benchmark studies.
Introduction
During the last decade great progress has been obtained in the understanding of the
locking of finite element methods for various "thin structures". For one-dimensional
problems, i.e. beams and arches, it is now completely known how "locking-free"
finite element methods should be constructed, cf. e.g. [1, 4]. For plates based on
the Reissner-Mindlin theory considerable advances have also been achieved, and by
now the problem with the shear locking can safely be claimed to be solved. In this
respect we refer to the papers [3, 8, 10, 48] in which the optimal order of convergence
is rigorously established for several methods and families of methods.
With regard to shells the present situation is unfortunately far from satisfactory.
Even if the literature on the field is voluminous, it is commonly admitted that the
elements presently in use are not completely reliable, cf. e.g. the survey [11]. From
the viewpoint of numerical analysis, it is not surprising to encounter difficulties with
shells. In fact, the general concept of "shell" covers whole families of problems with
very marked differences in behaviour depending on e.g. the geometry of the mid-
surface, the boundary conditions and the loading. Hence, the goal of developing
"the shell element" may be too ambitious, at least at present.
There are two main classes of shell behaviour that can be clearly distinguished:
the membrane and bending dominated cases. Mathematically, the membrane case
is a singularly perturbated second order elliptic problem when the thickness of the
shell is small. For this case the standard finite element method usually works quite
well. Difficulties arise in the other case, i.e. when the deformation is bending
dominated. Then, the limit problem leads to constraints which give rise to locking
if they are exactly imposed in the finite element model.
The numerical analysis of the locking problems for shells is still in its infancy.
The papers are in fact so few that they can very well be briefly reviewed here.
Pitk-aranta was the first to address this question in [40], where he considers hp
finite element methods for a cylindrical shell in a bending dominated state. He
shows that the standard h-version locks for low order methods, or if the finite
element mesh is not aligned with the axis of the cylinder. In the paper it is also
shown that the p-version with a fixed mesh is free from locking. Finally, a carefully
designed hp-method is shown to be asymptotically convergent both with respect
to h and p, but unfortunately only on an aligned rectangular mesh. Pitk-aranta
analyses the method directly from the (modified) energy expression. The "partial
selective reduced integration" introduced in the more recent paper by Arnold and
Brezzi [2] is based on the classical approach of writing the problem in mixed form.
They first use a splitting of the energy (already used by Pitk-aranta in [40]) to write
the mixed system in such a form that the "Z-ellipticity" condition of the classical
saddle point theory is avoided [9]. Then they prove the "inf-sup" condition by
constructing a Fortin operator. For this, the technique of "bubble functions" is
used. In the construction of the Fortin operator, the authors need the assumption
that the geometrical parameters (i.e. the fundamental forms and the Christoffel
are piecewise constants, which severely restricts the applicability of their
results (namely to circular cylinders). It seems to be non-trivial to extend their
analysis to the general case. In [49] Suri studies the approach of Arnold and Brezzi
as an hp-method for rectangular elements. The analysis of Suri is in the spirit of
Pitk-aranta, and the modified energy expression is analysed directly without the use
of the equivalent mixed method. Also here the assumption of piecewise constant
coefficients is used. Let us remark that, intuitively, this assumption seems to be
even more restrictive for a "p-dominated extension" procedure as the elements then
can be quite large and the variation in the geometry big. Let us finally mention the
paper by Kirmse [28] in which a spherical surface is considered and a locking free
method is designed.
The purpose of the present paper is to explore the application of stabilization
techniques to shell problems. These stabilization techniques have earlier been shown
to give methods free from locking for several related problems in continuum mechanics
such as incompressible elasticity [26, 22, 21], beams [31], arches [30] and
Reissner-Mindlin plates [39, 27, 47, 33]. The stabilization technique is quite simple.
The equations are first written in a variational mixed form. Then, properly weighted
least-squares-type expressions of the equilibrium and constitutive equations are
added to the bilinear form. This is, indeed, how we here treat the equations of the
Naghdi shell model [35, 5]. For the resulting finite element method we are able to
prove the stability, hence the method is free from locking. For our result we do not
need any restrictive assumptions on the geometry. Neither do we need any stabilizing
bubble degrees of freedom, and we are able to use standard finite element
spaces. we should also emphasize that the the auxiliary variables introduced in
the mixed stabilized bilinear form are condensed in an implementation which then
takes the standard displacement form. Our method is formulated in the h-version.
In principle, the method can be formulated in a general hp setting, but as we extensively
use inverse estimates the resulting method would not be uniformly stable
with respect to the polynomial degree p. This lack of stability might, however, be
compensated for by the better approximability properties of an hp-method, especially
if a properly refined mesh is used. A general analysis would require that all
technical results are verified also with respect to p. As the present paper is already
rather technical, we have not tried to perform this.
In other respects our analysis suffers from the same shortcomings as the previous
papers on the subject. In order to concentrate on the locking, we give our
estimates assuming that the exact solution is sufficiently smooth. As it is known
that the shell problems have boundary layers of different length scales [41], this
assumption is not very realistic. But as we have a stable method, the treatment of
the boundary layers should be done by a proper mesh refinement and this is a problem
of approximation theory. The biggest shortcoming for the method we present,
and also for the methods that have been mathematically studied in previous works
[2, 49, 28], is that they do not perform properly when they are applied to membrane
dominated shells. Hence, we are still far from having "the shell element".
An outline of the paper is as follows. In the next section we specify our notation
regarding the differential geometry of shell surfaces. In Section 3 we recall the
Naghdi shell model and give the equations in the bending dominated case that we
focus on. In the next two sections we introduce and analyse two finite element
formulations, corresponding to different choices of the stabilization weights. Then,
in Section 6, we give numerical results obtained with these methods in benchmark
computations. Finally, we close the paper with an appendix in which we prove some
results of differential geometry that we have used in our analysis.
Shell geometry and notation
we use the classical representation of the shell geometry, as described in [24]. Thus,
the shell mid-surface is characterized by a map ~r which is a one-to-one mapping
from
\Omega\Gamma an open domain of IR 2 , into IR 3 . ~r will be assumed to be as smooth as
required (in practice C 3 ( -
\Omega\Gamma22 The actual surface we denote by S, i.e.
-\Omega\Gamma4 we
consider a shell of uniform thickness t.
Regarding the concepts and quantities of the differential geometry of surfaces
that we use in this study, we again refer to [24] for both definitions and notation.
However, with a view to making equations more compact, we also make use of
the alternate notation by which surface tensors are represented by letters with a
number of underbars equal to their order. In particular, a scalar is denoted by
a simple letter. With this notation, a tensor is considered as independent of the
type (covariant or contravariant) of its representations in the curvilinear coordinate
system, and transforming one set of components into another is simply done by
using the transformation formulas. Thus we will write
ff\Omega ~a
ff\Omega ~a
ff ~a
ff\Omega ~a
\Deltafi ~a
Here ~a are the covariant surface base vectors and ~a ff
is the contravariant base. A single dot will denote the simple contraction operator:
ff\Omega ~a
ff\Omega ~a
whereas a colon denotes a double contraction:
The metric tensor for the surface we denote by a and the second and third
fundamental forms are denoted by b and c (= b : b), respectively.
For a second-order tensor, we define the transposition and symmetry operators
as follows
we next consider differential operators. The gradient of a tensor is obtained by
taking the covariant derivative and adding one index
and the divergence by contracting the gradient on the last two indices
For a surface integral over !, a subdomain of \Omega\Gamma we use the following compact
notation Z
Z
a d- 1
we also write Z
ds def
Z
ds
for an integral along a curve C, with (ds) . With this notation, Gauss
theorem reads Z
Z
where - is the unit outward normal, in the tangential plane, to the boundary ~r(@!).
From this the following Green formulas are derived
Z
Z
ds \Gamma
Z
(div u)w dS; (2.2)
Z
Z
Z
For our analysis it will be convenient to use special norms for Sobolev spaces
of tensors. If two zero-order tensors v and w are in L 2 (!) we define the "intrinsic"
inner-product and norm
Z
This new norm is equivalent to the regular L 2 (!)-norm, cf. [43]. we also extend
this definition to higher order tensors and write
Z
Z
we also denote the corresponding norms by k \Delta k 0;! and these can be shown to be
equivalent to the regular L 2 (!)-norm of any of the representations of these tensors in
covariant or contravariant components [43]. Finally, we introduce new H 1 (!)-norms
through
Recalling the relationship between covariant and regular derivatives, it is clear that
this new norm is again equivalent to the standard H 1 (!)-norm of the the tensors
in any representation. Using higher-order covariant differentiation, there is no
difficulty in similarly defining intrinsic H k -norms (k - 2) and establishing the corresponding
equivalence properties. For all norm and inner-product signs, we will
omit the domain subscript when the domain considered
is\Omega itself.
3 The variational shell model
The shell model we consider, the so-called Naghdi model, is of the Reissner-Mindlin
type, i.e. it includes the effect of shear deformation, cf. [35, 5, 50, 52, 19] and the
references therein. The unknowns are the displacement at each point
of the mid-surface decomposed into tangential and transverse parts (respectively a
first-order and zero-order surface tensor), and ' as the first-order tensor representing
the rotation of a fiber normal to the mid-surface S in the undeformed configuration.
we let
we assume that the boundary conditions on
are given so that ~r(\Gamma 0 ) is the part of the boundary along which the shell is
fully clamped (i.e. the part where all displacements are left
free. The problem is posed in the
domain\Omega and we define the displacement space
U by
The variational formulation then reads [5]: find (~u; ') 2 U such
\Omega
Z
\Omega
Z
\Omega
Z
\Omega
Here the material properties are given by the two tensors
a ff- a fi- + a ff- a fi-
a fffi a -
and -, ", fl denote the bending, membrane and shear strain tensors, respectively:
we also recall that b and c are the symmetric tensors corresponding to the second
and third fundamental forms of the surface, respectively.
we also write equation (3.1) in the shorter form:
with
Z
\Omega
and
Z
\Omega
Z
\Omega
Implicit in the use of a shell model is the assumption that the shell is "thin".
Hence, it is essential to study how the solution behaves in the limit when the
thickness This limit behavior is very different depending on the geometry
of the shell and the boundary conditions and depends in a crucial way (cf. [44, 45,
15, 16, 36, 37]) on the subspace
consisting of the displacement fields with vanishing shear and membrane strains.
In this paper we only consider the case when the shell is said to be in a bending
dominated state of deformation. This is the case when (cf. [44, 16, 37])
Then the limit problem is obtained by assuming the loading to be given by
f
with ~
f independent of the thickness t. In the sequel we therefore study the following
problem.
U such that
Here we introduced the notation F (\Delta) def
From [6] and [18] we know that, for any t ? 0, P t has a unique solution, provided
that the distributed force field ~
f is in an appropriate space, say [L
Furthermore, as t tends to zero, P t can be seen as a standard penalty problem. The
limit problem P 0 is defined as:
such that
The term "limit problem" is justified since, as t ! 0, it can be shown that
the solution (~u; ') of P t converges strongly in U to the solution (~u
[13, 40]).
In the approximation of plate problems it has often turned out to be advantageous
to use a mixed formulation, cf. [9, 10, 20]. This is the approach that we use
here as well. Hence, we define as new unknowns stress variables dual to the strains
in the variational formulation (3.5). By n we denote the symmetric membrane force
tensor, Here we depart from the usual notation, cf. Remark 3.1 below. The shear
force we denote by q. These are connected to the membrane and shear strains
through the constitutive equations
The symmetric bending moment tensor is not taken as an independent unknown,
but we use the following abbreviation based on the bending constitutive equation:
By introducing
as new unknowns and by writing the constitutive equations (3.8) and (3.9) in a
weak form, we obtain an equivalent "mixed" formulation of the problem:
Q such that
with the new bilinear forms defined through
Z
\Omega
and
Z
\Omega
Here we use the notation
a ff- a fi- + a ff- a fi- \Gamma 2-
a fffi a -
a
The second order tensor -
G is the inverse of G (i.e. G \Delta -
G
is such that, for any two symmetric tensors T and
Finally, we recall that using the symmetry of the second order tensors, and
the Green formulas (2.2) and (2.3), the following integration by parts formulas are
obtained
Z
Z
Z
div (b \Delta m) \Delta v dS (3.12)
Z
(b
Z
(m ds \Gamma
Z
Z
Z
Z
Z
(n ds
and
Z
Z
Z
(b
Z
div q v 3 dS+
Z
Collecting these formulas gives
Z
Z
Z
Z
Z
Z
Z
(m ds (3.15)
Z
ds
Z
Hence, recalling the constitutive equations (3.7)-(3.9), problem P t gives the
following equations of equilibrium
in the
domain\Omega\Gamma and the natural boundary conditions
on the
Remark 3.1 In the above equations the quantity
is the tensor obtained by integrating the three-dimensional membrane stresses over
the shell thickness, cf. e.g. [29, 24, 35, 52]. Usually this is called the membrane force.
For designing a stable method it seems more practical to introduce the variable n,
and not ~ n, as an new unknown. For simplicity we refer to n as the membrane force.
4 A first approximation scheme
we use a finite element partitioning C h of
-\Omega into straight-sided triangles or quadri-
laterals. Naturally, the partitioning is assumed to satisfy the usual regularity and
compatibility conditions, cf. [14, 7, 42]. The diameter of an element K 2 C h we
denote by hK , and we let
hK . By \Gamma h we denote the collection of edges
in the mesh, and by hE the length of an edge . The letters C and c are
henceforth used to denote generic strictly positive constants, independent of both
t and h, which are allowed to take different values at different occurrences except
when appearing with indices.
Let U h and Q h be finite element subspaces of U and Q, respectively. A mixed
finite element method of the "classical" saddle point type would be based on the
formulation M t , viz.
such that
It appears difficult to design a stable finite element method directly based on this
formulation. Hence, we follow an approach that has turned out to be fruitful for the
related beam, arch and plate problems, cf. [22, 21, 47]. The approximation scheme
we propose is derived from (4.1) by adding weighted least-squares-type terms of
the equilibrium equations (3.16)-(3.18) and the constitutive relations (3.7)-(3.8).
In order to simplify the presentation of our method, we introduce the weighted
inner-product
with analogous definitions for higher order tensors. Further, in accordance with
and (3.22) we introduce the notation
Our first stabilized method will now be defined as follows:
such that
where
with a stabilizing term originating from the equilibrium equations
div
\Gammaff 3
and one from the constitutive equations
Z
\Omega
dS (4.7)
Z
\Omega
dS:
The right hand side is defined as
The norm used for the displacement is the modified H
1(\Omega\Gamma9744/ that we previously
defined, i.e.
For the stresses, we use a discrete norm defined by
I
with "interior" terms
and "jump" terms
Z
ds
where stands for the jump of the quantity if the edge is in the interior, and for
the quantity itself if it is on \Gamma 1 . Finally, we let
t;h
we begin our analysis by noting the consistency of formulation S h
t with respect
to the continuous problem M t .
Lemma 4.1 Suppose that ~
2(\Omega\Gamma7 3 . Then the solution (~u; '; n; q) 2 U \Theta Q of
First we note that, due to the assumption ~
3 , the quantity
r) is well defined for all (~v; j; p; r) 2 U h \Theta Q h . Therefore, using
the continuous formulation M t , the equilibrium equations (3.16)-(3.18) and the
constitutive equations (3.7)-(3.8), we directly obtain
div
div ~
\Gammaff 3
Z
\Omega
dS
Z
\Omega
dS
In order to obtain a stable method we have to specify the finite element spaces.
To this end we let, for l - 0,
l (K) if K is a triangle,
and for k - 1 we now define
and r ff j K 2 R
Remark 4.1 With the choice of finite elements spaces that we made, the stresses
are interpolated discontinuously between elements. This implies that they can be
eliminated at the element level to obtain a purely displacement-based numerical
scheme, cf. equation (6.3) below.
Remark 4.2 Our choice of covariant components for the displacement variables,
and of contravariant ones for the stresses, enables us to establish Lemma 4.4 (using
Lemma A.2), which we need for proving the stability of the method.
we are now ready to state our main stability result. we will assume that the
shell thickness is in the range
Lemma 4.2 Assume that:
where the positive constants C I
are fixed, derived
from inverse inequalities.
(S4) h is sufficiently small.
Then the stabilized formulation S h
t is stable, i.e. there is a positive constant C such
that, for all ( ~
there exist (~z; j;
and
be arbitrary. The proof we will now divide into
three steps.
Step 1. First, a direct calculation gives
Z
\Omega
Z
\Omega
Z
\Omega
Z
\Omega
we will use the following inverse estimates
For clarity, we will establish (4.21) and (4.22) as a separate result in Lemma 4.3
below. we also have
and from Lemma A.1:
Z
\Omega
so that we get X
Combining equations (4.21), (4.22) and (4.25) with (4.20), we obtain
Z
\Omega
Z
\Omega
Z
\Omega
Z
\Omega
the ellipticity result proved in [6, 18] gives
Z
\Omega
Z
\Omega
when choosing
and for h small enough. Next, Lemma A.1 implies
Z
\Omega
Further, it holds Z
\Omega
Since
a combination of the above relations gives
again when h is small enough, to take care of the term \GammaC h 2 k~w; -k 2
0 in (4.26).
Finally, using the arithmetic-geometric mean inequality, and the boundedness of b,
we get, for
Hence, by setting j close enough to one, we obtain
I
Step 2. Next, we use the result that there exists (~v; such that
Z
ds
and X
The proof of this result is rather technical and hence we will postpone it to Lemma
4.4 below.
With this ~v we first write
with
I
Z
\Omega
Z
\Omega
div
div ~
\Gammaff 3
and
Z
\Omega
Z
\Omega
Note that estimates (4.21) and (4.22) immediately imply, for all (~x; -) 2 U
In order to bound I, we apply these relations to both ( ~
and we use
the Schwarz inequality and (4.36). we obtain
I
I
we transform the term II by using the integration by parts formulas (3.13) and
over each element
Z
Z
with
Z
ds
Z
Hence, using (4.31), (4.35), (4.36) and the arithmetic-geometric mean inequality we
get
Combining (4.37)-(4.44) and using, once again, the arithmetic-geometric mean in-
equality, we obtain
I
I
2oe
I
for any strictly positive oe. By choosing oe small enough, we then have
I
Step 3. Choose now (~z; j;
relations
(4.34) and (4.45) proved in the previous steps then give
I
when choosing Using (4.36) we get
and the assertion is thus proved.
Next, we will prove the missing steps in the proof above. we start with the
inverse estimates (4.21) and (4.22).
Lemma 4.3 The following inequalities hold for any (~v;
First we note that
while, by virtue of the smoothness of b we also have
Therefore the proof will be identical for both estimates. Further, the above definition
of Sobolev norms for tensors implies
where the norm in the right hand side now denotes the standard H 1 -norm. Our
argument will not require the detailed expression of m fffi (~v; j), so we instead symbolically
denote
where all coefficients involved are smooth functions that incorporate geometric and
material terms. For a smooth function ae we now define T h ae as the linear part
of its Taylor expansion at the center of each element. Using standard expansion
properties, we have for any (ff; fi) and any element
k(ae fffii-
k(ae fffii-
Hence, we obtain
we can now invoke standard inverse estimates on polynomial functions to obtain
and, using (4.53) combined with an inverse estimate, we finally get
which, combined with (4.24), concludes the proof.
Remark 4.3 It is not clear whether estimates (4.47) and (4.48) would hold without
the terms h 2 k~v; jk 2
0 . These estimates, however, are sufficient for our purposes in
their present form.
Remark 4.4 In the above proof, it clearly appears that estimates (4.47) and (4.48)
follow from similar local estimates that we can explicitly write as:
This implies that these inverse inequality constants can be computed for each element
using a simple local eigenvalue problem, cf. [25] where this idea is introduced
for several stabilized methods. Therefore the choice of ff 1 and ff 2 can be made
automatically in the course of the assembling process.
we proceed to establish the result required at the second step of the proof of
Lemma 4.2. For this, the technical Lemmas A.1 and A.2 are essential.
Lemma 4.4 There is a positive constant C such that, for all (k; s) 2 Q h , there
exists (~v; such that
Z
ds
ds - Cjk; sj J
and X
Proof : we use the normal Lagrange degrees of freedom for ~v and we choose these
so that at the vertices and the internal nodes of all elements. To the remaining
degrees of freedom of ~v along edges we assign values so that
where a bar over the symbol of a continuous function denotes its value at the
midpoint of the edge, and b E stands for the second-degree "bubble function" along
the edge with value one at the midpoint (i.e. - b As it holds
this construction yields a function ~v lying in the appropriate finite element subspace
for the displacement.
Now, we first note that Lemma A.2 gives
Z
Z
ds
Z
where i is a smooth positive function bounded away from zero. Then we can apply
Lemma A.1 to obtain
Z
ds - C
Z
ds
Z
and since the quantities involved are now polynomial functions, a standard scaling
argument yields
Z
ds - C
Z
ds
so that, using again Lemmas A.1 and A.2, we have
Z
ds - C
Z
ds
Z
ds
To bound v, we first note that inverse estimates imply
Then, using the fact that all degrees of freedom are equal to zero except along edges,
we get by scaling arguments, and Lemmas A.1 and A.2
For the component v 3 similar arguments give
Z
ds - C
and
The estimates (4.62)-(4.65) then prove the assertion.
To perform the error analysis, interpolation estimates in non-standard norms
will be required.
Lemma 4.5 For (~u;
[H k+1
it holds
Proof : The estimate follows from scaling arguments.
Lemma 4.6 For (n;
[H
it holds
r) be the L 2 -projection of (n; q) in Q h . Since no continuity is
imposed on the finite element functions, this projection is defined locally element
by element. we now have
The last term above is handled by a scaling argument
and the asserted estimate then follows from standard estimates.
Let us now prove the error estimate for our numerical scheme. In the proof we
will repeatedly need the estimate
which is a direct consequence of standard inverse inequalities.
Theorem 4.1 Assume that conditions (S1)-(S4) of Lemma 4.2 hold and that the
solution of M t is such that ~u 2 [H
and q 2 [H k
. Then the finite element formulation S h
t has a unique solution
satisfying
be the interpolant to (~u; '; n; q) satisfying
the estimates of Lemmas 4.5 and 4.6. By the stability result of Lemma 4.2 there
exists (~v; j; p; r) 2 U h \Theta Q h such that
and
Next, the consistency of Lemma 4.1 gives
To proceed, we write out the expression for the bilinear form
For the first term above the Schwarz inequality, Lemma 4.5 and (4.70) give
By using the integration by parts formulas (3.13), (3.14) and Lemma 4.5 we get
Next, we directly get
and
The first stabilizing term is estimated using Lemmas 4.5 and 4.6 as follows
\Theta
where we used the inverse inequality (4.68). For the second stabilizing term we
obtain
by virtue of the previous estimates (4.75)-(4.77), the Schwarz inequality and Lemma
4.5. Finally, we note that Lemmas 4.5 and 4.6 also give
Hence, collecting (4.71) to (4.80) and using the triangle inequality we get the asserted
estimate.
we will close this section by giving some remarks on the analysis of the method.
In [2] the analysis is performed using an abstract semi-norm for the stresses, instead
of the discrete norm we have used. This abstract semi-norm is defined as follows
(~v;j)2U
For our finite element method it is possible to prove the stability using the same
global norm as in [2], i.e.
This is achieved using the "Pitk-aranta-Verf?rth trick" introduced by Pitk-aranta in
connection with Babu-ska's method for approximating Dirichlet boundary conditions
[38] and Verf?rth for mixed methods for the Stokes equation [51]. The idea has later
been extensively used for both classical and stabilized formulations of the Stokes
problem, cf. [46, 21].
we are, however, not able to carry through the whole error analysis with this
norm. The reason is that the stabilizing term (4.78) which we are not able to bound.
One possibility for an analysis would be to use the combined norm jjj \Delta jjj +k \Delta k t;h
for the stresses. This is possible since the same stability construction can be shown
to yield the stability with respect to both norms for the stresses. we choose a
different approach and show that the estimate for the stresses in the norm (4.81)
can be proven "a posteriori".
Theorem 4.2 Under the assumptions of Theorem 4.1 it also holds
From the definition of the semi-norm (4.81) there is (~v;
that
and
Now we let (~v be the Cl'ement interpolant to (~v; j). From [17] we have
Next, we write
For the first term above an integration by part using (3.13) and (3.14) gives
From the consistency we have
Hence, it holds
\GammaS c
For the first term above we get
The second term is estimated as follows
\Theta
\Theta k~v h
where we used the inverse inequality (4.68). For the third term we directly get
Combining the above estimates gives
From Lemma 4.5, the inverse estimate (4.68) and Theorem 4.1 it follows that
Hence, the assertion follows from Theorem 4.1.
5 A second stabilized method
In this section, we propose an alternative stabilized scheme, obtained from the
previous formulation by changing the weights in front of the stabilizing terms that
originate from the constitutive equations, in the spirit of what is considered for
plates in [39, 27, 47, 34, 33, 32]. This second stabilized method reads:
~
such that
~
where
~
with the stabilizing term S e
h originating from the equilibrium equations defined as
before in (4.6), and the term from the constitutive equations redefined by
~
Z
dS
Z
dS:
The right hand side F h remains unchanged as defined by (4.8).
For the displacement variables we will use the following norm
which is obviously stronger than the H 1 -norm used for our first method. The norm
for the stresses we now have to weaken by first redefining the "jump" term
Z
ds
and then defining
I
by keeping the interior part as defined in (4.11). we will also use the combined
notation
t;h
The analysis of this new scheme will now essentially follow the same lines as in
Section 4, so there is no need to give the full proofs for all results. we will instead
highlight the differences between the two formulations, whenever these differences
are significant.
First, we again note that the formulation is consistent.
Lemma 5.1 Suppose that ~
2(\Omega\Gamma7 3 . Then the solution (~u; '; n; q) 2 U \Theta Q of
~
The finite element subspace U h \Theta Q h we keep as defined by (4.16) and (4.17) in
Section 4. With this choice we can state the new stability result.
Lemma 5.2 Assume that:
where the positive constants C I
are fixed, derived
from inverse inequalities.
sufficiently small.
Then the stabilized formulation ~
t is stable, i.e. there is a positive constant C such
that, for all ( ~
there exist (~z; j;
~
and
Proof : we again go through 3 steps.
Step 1. First, a direct substitution yields
~
Z
w)dS
Z
Note that, for any K,maxf1; ff i
Using as before estimates (4.21), (4.22) and (4.25), we thus obtain
~
we again set
C I
2:
Then, proceeding like in the proof of Lemma 4.2, we can conclude that
~
I
Step 2. By changing the weights from hE to h 3
in (4.60), Lemma 4.4 changes as
follows: there exists (~v; such that
Z
and X
we first note that
By substituting this ~v we get
~
with
I
Z
Z
div
div ~
\Gammaff 3
Z
\Omega
Z
\Omega
and
Z
Z
Using the Schwarz inequality and estimates (4.40)-(4.41), then (5.12), we obtain
I
I
Further, recalling (5.7), we have
Z
Z
s
Z
s
Z
using the Schwarz inequality and (5.12). In the second term we again integrate by
parts using (3.13) and (3.14) , which gives
with
Z
and
Z
Z
Equation (5.11) implies
Hence, recalling (5.10), the arithmetic-geometric mean inequality gives
As in the proof of Lemma 4.2, we now combine (5.13)-(5.22) and use the arithmetic-geometric
mean inequality to obtain
~
I
Step 3. The assertion again follows from (5.9) and (5.23) by choosing (~z; j;
we now examine the modifications in the interpolation estimates. In view of the
stronger norm used for the displacements, a stronger assumption on the regularity
of ~u is required.
Lemma 5.3 For (~u;
[H k+2
it holds
Using the definition of the norm [] \Delta [] t;h , we have
Therefore, recalling that the interpolation space is one degree higher for ~u than for
scaling arguments imply
The next term in (5.24) is also treated by standard scaling arguments, and the
remaining terms are unchanged from Lemma 4.5.
By contrast, the new norm for the stresses is weaker than the previous one, so
the interpolation estimates remain valid.
Lemma 5.4 For (n;
[H
it holds
we can now state and prove our final approximation result.
Theorem 5.1 Assume that conditions ( ~ S1)-( ~ S3) of Lemma 5.2 hold and that the
solution of M t is such that ~u 2 [H
and q 2 [H k
. Then the finite element formulation ~
t has a unique solution
satisfying
be the interpolant of (~u; '; n; q) satisfying
the estimates of Lemmas 5.3 and 5.4. By the stability result of Lemma 5.2 there
exists (~v; j; p; r) 2 U h \Theta Q h such that
and
Next, the consistency of Lemma 5.1 gives
~
To proceed, we write out the expression for the bilinear form
~
For the first term above, the Schwarz inequality, Lemma 5.3 and (5.27) give
By using the integration by parts formulas (3.13), (3.14) and Lemma 5.3 we get
Next, we directly get
and
The first stabilizing term is estimated using Lemmas 5.3 and 5.4 as follows
\Theta
where we used the inverse inequality (4.68). For the second stabilizing term we
obtain
by virtue of the previous estimates (5.32)-(5.34), the Schwarz inequality and
Lemma 5.3. Finally, we note that Lemmas 5.3 and 5.4 also give
Hence, collecting (5.28) to (5.37) and using the triangle inequality we get the asserted
estimate.
6 Numerical results
6.1 Numerical procedure
we have implemented the two stabilized methods analysed in Sections 4 and 5, for
the lowest-order finite element spaces, i.e. letting in (4.16) and (4.17), in the
simplified case of circular cylinders, using the MODULEF library. In addition to
the obvious implementational simplification, circular cylinders indeed feature two
major advantages:
ffl They provide one of the very few instances of shell geometries for which some
reference solutions can be derived, either in closed form or with arbitrary
numerical precision, so that reliable benchmarks are available.
ffl This geometry allows a wide variety of possible limit behaviours when the
thickness is very small, according to how boundary conditions are imposed. In
particular, if essential boundary conditions are not imposed outside of rulings
(straight lines parallel to the axis), it is easily seen that the limit problem is
bending-dominated.
we considered two different benchmarks. In the sequel, we refer to the natural
coordinate system, shown in Figure 1, in which ~a 1 and ~a 2 are unit vectors, respectively
tangent and normal to generators. The benchmarks are defined as follows:

Figure

1: The natural coordinate system
a) The boundary conditions and the loading are defined in such a way that the
shell behaves like an arch. Considering a cylindrical "slice" at the ends of
which one imposes assuming that the loading ~
f is a function
of - 2 such that f indeed easily seen that the resulting behaviour is
that of an arch, i.e.: (
we used an example, borrowed from [12], for which a closed-form arch solution
can be derived. This example, described in Figure 2, corresponds to a semi-circular
arch clamped at both ends and loaded by a uniformly-distributed
constant force. we considered a slice of length equal to 1.5 times the radius
R and, for symmetry reasons, one half of this slice was computed, so that the
actual computational domain was a 1:5 \Theta -rectangle.
b) A fully-circular cylinder is loaded by a periodic pressure, so that the shell
reduces to a one-dimensional model along the axis. This problem is analysed
in detail in [41], and we applied the procedure described therein to obtain
numerical solutions with arbitrary precision, using symbolic calculus software.
The specific example considered is a cylinder with free ends and of length 2R,
loaded by a pressure Figure 3. Due to symmetry,
only one sixteenth of the structure was effectively computed, using a 1 \Theta -rectangle.
In each case, numerical solutions were computed using triangular meshes automatically
generated by a Vorono-i method available in MODULEF [23], for gradually
~
computational
domain
R
1:5R
clamped

Figure

2: Semi-circular arch
R
computational
domain

Figure

3: Cylinder loaded by periodic pressure
MODULEF chapelle
sinus.nopo
547 POINTS
COIN BAS GAUCHE
COIN HAUT DROIT

Figure

4: Mesh generated for the 1 \Theta -domain (N=20)
refined discretizations of the rectangular boundary, dividing each side into N equal
parts, with One of these meshes is shown in Figure 4. we deliberately
avoided using meshes aligned with the sides of the domain, or with the axis of
the cylinder, that seem to alleviate locking phenomena in some instances [40].
The Poisson ratio - was set to 0.3. The stabilization constants were determined
by performing some preliminary calculations on the first benchmark. For the first
stabilized method, we then used ff
. For the second method we chose ff
In both cases the choice ff allowed because, with
each element for any discrete membrane stress field n (since all n fffi
are constant), the corresponding least-square term provides no further stabilization
than the first one.
In order to compute the stabilized numerical solutions, we used standard solvers
for symmetric positive matrices. To that purpose we eliminated the stress degrees
of freedom element by element, as indicated in Remark 4.1 above. Denoting by
U the column vector relative to the displacements and rotations, and by P that
relative to the stresses, the stabilized mixed methods indeed lead to the following
typical matrix equation:
PU
!/
U
FU
From the proof of the stability results, it is clear that both MUU and M PP are
symmetric positive definite matrices so that we can eliminate P from this equation.
Note that, since the stress finite element functions are discontinuous, this elimination
can be carried out at the element level as the assembling is performed. Thus,
instead of (6.2), the problem solved in practice is
where the matrix to be inverted, (MUU +M T
PP M PU ), is now symmetric positive
definite. Hence we solve a matrix equation similar to what would be obtained from
a standard displacement-based finite element method.
6.2 Analysis of the results

Figures

5 and 7 show relative errors in the H 1 -semi-norm for the displacements, i.e.
computed from solutions obtained by the first stabilized method, and also by a
purely displacement-based scheme corresponding to the standard Galerkin approximation
of problem P t using U h . Thickness values of 0.1, 0.01, 0.001 (scaled by R)
were considered. Figures 6 and 8 compare similar errors for the rotations with interpolation
errors, since optimal convergence is expected from the theoretical analysis.
These interpolation errors turned out to vary very little with the thickness (typi-
cally within 10% ranges), so for the sake of legibility only their values for are
displayed. Note that, from the theory, optimal error estimates are not expected for
the displacements since convergence is governed by the approximability properties
of the lowest degree finite element space, i.e. that of the linear rotations. This is
confirmed by the numerical results, as can be seen in Figures 5 and 7 where the
slope of the best convergence curves barely exceeds unity. This explains why we do
not plot interpolation estimates for the displacements.
In

Figures

5 and 6, i.e. for the arch benchmark, we see no significant influence
of the thickness on approximation errors of the stabilized method. Moreover, the
errors for the rotations are very close to interpolation errors. By contrast the
displacement-based method exhibits strong locking, with completely erroneous (in
fact vanishing) solutions for meshes except the finest one.
In the second example, a marked deterioration of convergence appears in Figures
7 and 8 for the first stabilized method when the thickness decreases. Approximation
errors are visibly affected in absolute value for each mesh, as well as in rates of
convergence. However, even errors for especially when
compared with the displacement-based method which dramatically fails here again.
we next compare the two stabilized methods. Figures 9 and 11 display the
relative errors for the displacements approximated by both methods for thickness
Stabilized
Stabilized
Stabilized

Figure

5: Benchmark a - Relative errors for ~u10.110 100
Stabilized
Stabilized
Stabilized

Figure

Benchmark a - Relative errors for '
Stabilized
Stabilized
Stabilized

Figure

7: Benchmark b - Relative errors for ~u10.110 100
Stabilized
Stabilized
Stabilized

Figure

8: Benchmark b - Relative errors for '
Stabilized 1
Stabilized 2
Reference t=0.1

Figure

9: Benchmark a - Compared relative errors for ~u10.110 100
Stabilized 1
Stabilized 2
Interpolation

Figure

10: Benchmark a - Compared relative errors for '
compared with the best numerical solutions obtained for thickness 0.1. Similar
errors for the rotations are plotted in Figures 10 and 12, with the corresponding
interpolation errors.
In the case of the arch benchmark, rotation errors remain essentially unchanged
in

Figure

10, but displacement errors now feature a more sensitive behaviour with
Stabilized 1
Stabilized 2
Reference t=0.1

Figure

Compared relative errors for ~u10.1
Stabilized 1
Stabilized 2
Interpolation

Figure

Compared relative errors for '
respect to the thickness parameter (Figure 9). However, this phenomenon is much
more limited than the deterioration observed with the first stabilized method in the
second benchmark. Furthermore, the numerical results obtained with this second
stabilized method here still strikingly differ from the locking behaviour exhibited
by the displacement-based method.
Regarding the results obtained for the second example and displayed in Figures
11 and 12, we see that the second stabilized method provides significantly improved
approximations for displacements as well as rotations.
In partial conclusion to our numerical tests, it appears that, even though both
stabilized methods give reasonably good results, the second one seems more promising
in that it limits the sensitivity of convergence behaviours with respect to the
thickness parameter in the two examples that we considered. This sensitivity is not
completely ruled out, but it does not compare with a real locking behaviour such
as the one displayed by the displacement-based method.
A

Appendix


This section contains results of essentially geometric character that are used in the
numerical analysis of the shell formulation. we start with a classical property of the
first fundamental form (see e.g. [43]), the proof of which we give for completeness.
Lemma A.1 There exist two strictly positive constants c and C such that, at any
point of -
\Omega and for any surface tensor j:
c
Consider the following function defined on R 2 \Theta -
(j
7! a fffi (M)j ff
Since a is the metric tensor of the surface, this function is strictly positive everywhere
except when (j it is continuous since ~r is smooth. Define
c and C as its minimum and maximum values on the compact set
Hence c ? 0 and (A.1) holds. Then (A.2) immediately
follows from the fact that the matrix (a fffi ) is the inverse of (a fffi ), so that
they have inverse eigenvalues.
we next establish a property regarding vectors orthogonal to curves which are
images in S of straight lines in \Omega\Gamma This result enters as a crucial element in the
stability of the proposed approximation schemes.
Lemma A.2 Let C be an oriented straight segment
in\Omega , - the unit tangent vector
of ~r(C) in S, - the unit normal vector (- any two
points on C, there exists a positive number i such that:
Moreover, i lies in an interval [i I ; i S ], with i I ? 0, the bounds of which are independent
of the specific segment considered.
parametrization of C. Define:
dx
By the chain rule:
and of course - equals ~ t up to a normalizing factor. Now let:
so that ~ t \Delta positively oriented. Therefore - equals ~n up to a
normalizing factor also, i.e.:
be any two points on C. we infer:
so that (A.3) holds. Moreover, from Lemma A.1 we have for any point M in \Omega\GammaC
c
Therefore:
r c
s
c
and setting i I =
C and
c , this completes the proof.
we now establish an ellipticity property for tensor -
Lemma A.3 For any symmetric tensor X at any point M
in\Omega , we have:
From the definition of -
in section 4, we have:
we can always construct, at least in a neighbourhood of M , a new coordinate
system such that, at M , the covariant base vectors are orthogonal and of unit-
length. we denote by ( ~
the components of X in this new coordinate system
where no difference subsists between covariant and contravariant forms. Since tensor
invariants are independent of the coordinate system considered, we have:
~
~
Thus, from (A.6):
since - 1
Acknowledgement

Most of this work was performed while R. Stenberg was
invited at INRIA by Professor M. Bernadou, under whose supervision D. Chapelle
was pursuing his doctorate. Both authors are grateful to Professor Bernadou for
giving them this opportunity, for several helpful discussions, and for his continuous
support during this project.



--R

Discretization by finite elements of a model parameter dependent problem.
Locking free finite elements for shells.
A uniformly accurate finite element method for the Reissner-Mindlin plate model
Mixed finite element methods for elastic rods of arbitrary geometry.
M'ethodes d'El'ements Finis pour les Probl'emes de Coques Minces.
Sur l'ellipticit'e du mod'ele lin'eaire de coques de
Finite Elemente.

Mixed and Hybrid Finite Element Methods.
analysis of mixed-interpolated elements for Reissner-Mindlin plates
Finite element analysis of shell structures.
A locking-free approximation of curved rods by straight beam elements
On the locking phenomenon for a class of elliptic problems.
The Finite Element Method for Elliptic Problems.
Analyse asymptotique des coques lin'eairement 'elas- tiques. I. Coques <Title>"membranaires"
Analyse asymptotique des coques lin'eaire
Approximation by finite element functions using local regulariza- tion

Mod'elisation des Coques Minces Elastiques.
A new finite element scheme for bending plates.
analysis of some Galerkin-least-squares methods for the elasticity equations
A new family of stable elements for nearly incompressible elasticity based on a mixed Petrov- Galerkin finite element formulation
MODULEF, Guide n
Theoretical Elasticity.
What is C and Inequalities for the analysis and design of finite element methods.
A new finite element formulation for computational fluid dynamics: VII.
A mixed finite element formulation for Reissner-Mindlin plate theory: uniform convergence of all higher-order spaces

A consistent first approximation in the general theory of thin elastic shells.
Stability, convergence and accuracy of a new finite element method for the circular arch pro- blem
Mixed Petrov- Galerkin methods for the Timoshenko beam
In preparation.
Stabilized MITC plate bending elements.
A stable bilinear element for Reissner- Mindlin plates
Foundations of elastic shell theory.
Energy estimates relating different linear elastic models of a thin cylindrical shell.
Energy estimates relating different linear elastic models of a thin cylindrical shell.
Boundary subspaces for the finite element method with Lagrange multipliers.
Analysis of some low-order finite element schemes for Mindlin- Reissner and Kirchhoff plates
The problem of membrane locking in finite element analysis of cylindrical shells.
Shell deformation states and the finite element method: a benchmark study of cylindrical shells.
Introduction 'a l'Analyse Num'erique des Equations aux D'eriv'ees Partielles.
Equilibre des Coques Elastiques Minces Inhomog
Statique et dynamique des coques minces.
Statique et dynamique des coques minces.
analysis of some finite element methods for the Stokes problem.
A new finite element formulation for the plate bending problem.
An hp error analysis of MITC plate elements.
A reduced constraint hp finite element method for shell problems.
The Nonlinear Theory of Shells through Variational Principles.

Mechanics of Solids with Applications to Thin Bodies.
--TR

--CTR
Douglas N. Arnold , Franco Brezzi , L. Donatella Marini, A family of discontinuous Galerkin finite elements for the Reissner-Mindlin plate, Journal of Scientific Computing, v.22-23 n.1, p.25-45, June 2005
Douglas N. Arnold , Franco Brezzi , L. Donatella Marini, A Family of Discontinuous Galerkin Finite Elements for the Reissner--Mindlin Plate, Journal of Scientific Computing, v.22-23 n.1-3, p.25-45, January 2005
Jaroslav Mackerle, Finite- and boundary-element linear and nonlinear analyses of shells and shell-like structures A bibliography (1999-2001), Finite Elements in Analysis and Design, v.38 n.8, p.765-782, June 2002
stabilized methods;naghdi model;shells;finite element methods