A Hybridizable Discontinuous Galerkin Formulation for
Non-Linear Elasticity
Hardik Kabariaa , Adrian J. Lewa , Bernardo Cockburnb
a Department
of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA, email:
[email protected] and [email protected]
b School of Mathematics, University of Minnesota, Minneapolis, MN 55455,email:
[email protected]
Abstract
We revisit the hybridizable discontinuous Galerkin method for non-linear elasticity introduced by S.-C. Soon [Ph.D. Thesis, University of Minnesota, 2008]. We show that
it can be recast as a minimization problem of a non-linear functional over a space of
discontinuous approximations to the displacement. The functional can be written as the
sum over the elements of the classic potential energy plus a new energy associated to
the inter-element jumps of the displacement. We then show that if this new energy is
not properly weighted, the minimizers might not converge to the exact solution. We
construct an example illustrating this phenomenon and show how to overcome it by suitably increasing the weight of the energy of the inter-element jumps. Finally, we explore
the performance of the method for the case of piecewise-linear approximations in rather
demanding situations in both two-dimensional and, for the first time, three-dimensional
situations. They include almost incompressible materials, large deformations with largeshear layers, and cavitation. We also compare the method with the continuous Galerkin
method and a previously explored discontinuous Galerkin method, and show that, when
using piecewise-linear approximations and a moderate number of degrees of freedom, the
current method turns out to be more efficient for the computation of the gradient.
Key words: Hybridized Discontinuous Galerkin, Non-Linear Elasticity.
1. Introduction
In this paper, we revisit the hybridizable discontinuous Galerkin (HDG) method for
non-linear elasticity introduced by S.-C. Soon [1, 2] and show how to rewrite it as a
minimization problem over finite dimensional spaces of discontinuous approximate displacements. We then show that, due to the discontinuous nature of the approximations,
the convergence of the minimizers to the exact solution can fail to take place in the general case. We then propose a simple way to avoid this phenomenon. Finally, we carry out
several numerical experiments displaying the performance of the method in several difficult situations in both two-dimensional and, for the first time, three-dimensional cases,
and with multiple different strain energy densities. We also evaluate the computational
efficiency of the resulting HDG method with respect to more traditional approaches, and
find very encouraging results.
Preprint submitted to Name of Journal
June 26, 2014
To better give a sense of the relevance of our results, let us place them in the context
of the development of discontinuous Galerkin (DG) methods for elasticity. In the last
decade, DG methods for continuum mechanics have been intensely explored. For example, DG methods have been proposed for linear elasticity in 2004 [3], for incompressible
and near incompressible elastic materials in 2002 [4] and in 2006 [5], for Timoshenko
beams in 2003 [6] and in 2006 [7], for fourth-order problems in structural and continuum
mechanics in 2002 [8] and 2005 [9, 10], for linear shells in 2007 [11], and for non-linear
elasticity in 2006 [12, 13]; see also the 2003 special issues in [14] and [15], and the
combination of XFEM and DG methods for the approximations of problems in fracture
mechanics proposed in 2008 [16]. An optimally convergent DG-method for XFEM was
later introduced in [17, 18]. DG methods for plasticity were considered in 2007 [19, 20].
The 2007 paper [21] and the 2008 paper [22] introduce a strategy to adaptively select the
stabilization parameters in non-linear elasticity. In 2007 [23], alternative ways of defining DG methods for non-linear elasticity are considered, in 2008 [24] a new DG method
for Kirchhoff-Love shells is studied, and in 2008 [25] a DG method for non-linear solid
dynamics is considered.
By the end of the 00’s, the HDG methods began to appear in this area. These methods
were devised so that they are amenable to static condensation. Thus, the structure
of the resulting stiffness matrices in these methods is somewhat special since, in our
setting, the only globally-coupled degrees of freedom are those of the approximation of
the displacement on the skeleton of the triangulation. The remaining degrees of freedom
are fairly inexpensive to solve for, basically involving static condensation, or the solution
of a small matrix involving the degrees of freedom associated with an element, for each
element in the mesh.
The first HDG method for linear elasticity appeared in 2009 [2]; its extension to nonlinear elasticity was proposed in 2008 [1]. In this method, the gradient of deformation and
the stress are computed more accurately than for other DG methods. Indeed, typically,
all previously mentioned DG methods use piecewise polynomial approximations of degree
k ≥ 0 for both the stress and the displacement, but, although the displacement converges
optimally, the stresses converge sub-optimally with order k. In contrast, the HDG method
[2] was experimentally shown to provide both approximate displacements and stresses
converging with the optimal order of k + 1, for k = 0, 1, 2, 3, as well as superconvergent
displacements with order k+2 when k = 2. Near the incompressibility limit, the orders of
convergence of the stresses were shown to degrade to 1.5 for k = 1, and to 2.75 for k = 2;
no superconvergence of the displacement took place for k = 2. Recently, the first error
analysis for this HDG method was carried out in [26]. For general polyhedral elements,
the orders of convergence of k + 1 for the displacement and k + 1/2 for the stress were
proven; in addition, the method was shown to be free of volumetric locking. The new
numerical experiments, carried out for triangular elements and uniform meshes, showed
that the orders of convergence found in [2] are actually not sharp. Indeed, for some
meshes and exact solutions, the orders of convergence for piecewise linear approximations
where experimentally found to be 1.5 for the symmetric part of the gradient, 1 for the
antisymmetric part of the gradient, and 2 for the displacement. Moreover, for quadratic
polynomial approximations, the displacement did not superconverge. All of these results
are consistent with those obtained in [2] for the near incompressible case. Only for
piecewise cubic approximations, the optimal convergence order of k + 1 = 4 for the
gradient and the displacement, as well as the superconvergence of order k + 2 = 5 for the
2
displacement did take place. These results still remain to be theoretically proven.
A second HDG method for linear elastodynamics was proposed in 2011 [27]; see also
[28]. The main difference with the first HDG method [2, 1] can be easily seen in the
case of linear, isotropic elasticity: σij = µ (ui,j + uj,i ) + λ uk,k . In contrast with the
first method, σ is not approximated. Instead, taking advantage of the identity σij,j =
µ ui,jj + pi where p := (µ + λ) uk,k , approximations of the gradient of the displacement
and of the hydrostatic pressure are provided. The resulting method converges with the
optimal order of k + 1 for both the displacement and the stress for k ≥ 0. Moreover,
superconvergence of parts of the approximate displacement allows for a postprocessing
resulting in a new approximate displacement of order k + 2 for k ≥ 1. These results were
proven in [29]. On the other hand, when the normal stress is imposed as a boundary
condition, the method in [28] loses the above-mentioned superconvergence, as shown in
[30].
A third HDG method for large deformation was proposed in [28] which does not reduce
to the HDG method considered in [27] (and [28]) in the linear case. In fact, such method
seems to be variation of the first HDG method for nonlinear elasticity introduced in [1],
as we show here. The numerical experiments presented in [28] for a neo-Hookean material
and piecewise polynomial approximations of degree k = 1, 2, 3, 4, indicate convergence
of order k + 1 for the displacement, the gradient of the displacement and the stress,
except for k = 1, 2 case in which the gradient converges with the suboptimal order of
k. Accordingly, the superconvergence of the displacement does not occur for k = 1, 2.
However, it does occur for k = 3, 4. These results are in full agreement with those
recently obtained for linear elasticity in [26]. This confirms the close relation between
the method proposed in [28] and in [2].
Finally, a fourth HDG method for linear elasticity was proposed in 2013 [31]. It was
constructed so that superconvergence of the approximate displacement holds, uses the
infinitesimal rotation as an additional unknown, and provides a weakly symmetric stress.
It has bigger local spaces for the stresses, and, more importantly, assumes that inverting the constitutive tensor is easy, which is not a reasonable assumption for nonlinear
elasticity as it requires the inversion of the first Piola-Kirchhoff stress tensor.
In this paper, we continue the ongoing effort on HDG methods for non-linear elasticity
and revisit the HDG method proposed in [2, 1] to obtain three new results. The first
is to show that the method can also be defined as a minimization problem of a nonlinear functional over a space of discontinuous approximations to the displacement, as
originally done in [12] for a DG method. Advantages of this type of formulations are that
the resulting stiffness matrices are symmetric, that numerical minimization strategies
could be used to solve for the approximations, that the variational structure mimic that
of the exact non-linear elasticity problem, and that this structure could be convenient
for the analysis of the method [32]. The functional to be minimized can be written
as the sum over the elements of the classic elastic potential energy plus a new energy
associated to the inter-element jumps of the displacement. The so-called stabilization
function controls the relative importance of these two energies in the functional, and
its proper choice is essential for the convergence properties of the method. Indeed, our
second contribution is to show that if the stabilization function is not properly chosen,
the minimizers might not converge to the exact solution at all. Indeed, we construct
an example illustrating this phenomenon and show that it is overcome by increasing the
weight of the energy of the inter-element jumps. A method to estimate the minimum
3
value of the weight of the inter-element jumps in the context of a non-hybridized DG
method for non-linear elasticity was introduced in [21, 22]; a similar or improved strategy
for HDG is yet to be devised. Our third contribution is of experimental nature and is
twofold. First, we explore the performance of the method for the case of piecewise-linear
approximations in rather demanding situations in both two-dimensional and, for the first
time, three-dimensional situations. They include almost incompressible materials, large
deformations with large-shear layers, and cavitation. We find that the method provides
very good approximations in all these cases. Second, we compare the method with the
continuous Galerkin (CG) method and a previously explored DG method [21, 22], and
show that, in our implementation for piecewise-linear approximations, the method turns
out to be more efficient for the computation of the gradient beyond a moderate number of
degrees of freedom. Comparisons in terms of the number of degrees of freedom, number
of non-zero entries in the stiffness matrix, and computing time are presented.
The organization of the paper is as follows. After formulating the non-linear elasticity problem as a minimization problem in §2, in §3 we define the HDG method as a
minimization problem and show that it is equivalent to the method proposed in [2, 1]. In
§4, we show how, when the stabilization function is taken too small, convergence of the
minimizers to the wrong solution can actually take place. In §5 we present the numerical
results, and end in §6 with some concluding remarks.
2. The Non-Linear Elasticity Problem
We consider a body with reference configuration B0 ∈ R3 that deforms under the
action of external loads. We assume that B0 is an open, connected Lipschitz domain.
Due to the deformation, a point X ∈ B0 is mapped to a point x = ϕ(X) in the deformed
configuration, where ϕ : B0 → R3 is the deformation mapping. Here we are dealing
with simple linear and non-linear elastic bodies, namely, bodies made out of hyperelastic
materials for which there exists a strain-energy density function W : B0 × R3×3
→ R,
+
W (X, F), where F = ∇ϕ(X) is the deformation gradient at point any X ∈ B0 and
R3×3
is the set of 3 × 3 matrices with positive determinant. For such materials, the
+
first Piola-Kirchhoff stress tensor is defined as P = ∂W/∂F. The strain-energy density
function W is assumed to satisfy the postulate of material frame indifference [33], which
states that the strain-energy density function should be invariant under superimposed
c (X, C), where C = FT F is the
rigid body rotations. Hence, there exits a function W
c (X, C).
right Cauchy-Green strain tensor, such that W (X, F) = W
The body is subjected to imposed deformations on ∂d B0 ⊆ ∂B0 , so that
ϕ = ϕd
on ∂d B0 ,
(1)
for some prescribed ϕd : ∂d B0 → R3 . For simplicity, we assume that ∂d B0 is a relatively
closed set with non-empty interior in ∂B0 . Static equilibrium configurations of the body
are rendered by stationary points of the potential energy functional
Z
Z
I[ϕ] =
W (X, ∇ϕ) − f · ϕ dV −
T · ϕ dS.
(2)
B0
∂N B 0
among admissible deformation mappings, namely, those ϕ ∈ V = {ϕ ∈ W 1,1 (B0 , R3 ) |
ϕ = ϕd on ∂d B0 , det ∇ϕ > 0 a.e. in B0 }. Here f : B0 → R3 is the body force per unit
4
volume, and T : ∂N B0 → R3 are tractions applied on ∂N B0 = ∂B0 \ ∂d B0 . Both f and T
are assumed smooth enough so that the integrals in (2) are well defined for any admissible
deformation mapping. The weak form of the Euler-Lagrange equations is
∂
0 = hδI[ϕ], δϕi =
I[ϕ ]
∂
=0
Z
Z
(3)
=
T · δϕ dS,
P(X, ∇ϕ) : ∇δϕ − f · δϕ dV −
∂N B0
B0
where ϕ = ϕ + δϕ, for all smooth admissible variations δϕ such that δϕ = 0 on ∂d B0 .
Here A : B = trace(AT B), for any A, B ∈ R3×3 .
Only local minimizers of the potential energy functional described in (2) are stable
under small perturbations, and therefore are often the solutions of interests. Here we
assume that conditions for minimizers to exists are met, such as the strain-energy density
function W being polyconvex (see [34]). In particular, we will consider strain-energy
density functions of the Mooney-Rivlin type extended to the compressible range,
W (F) = U (det F) +
µ1
µ2
F: F +
(det F)2 F−1 : F−1 ,
2
2
(4)
where
U (J) =
λ
(log(J))2 − µ1 log(J),
2
(5)
and λ, µ1 , and µ2 are material constants. Since we shall only consider homogenous
bodies, we have dropped here the dependence of W on X. When µ2 = 0, the strainenergy density function for the neo-Hookean material extended to the compressible range
is recovered. For completeness, the expression for the first Piola-Kirchhoff stress tensor
obtained from (4) is
P(F) = λ log(det F) − µ1 + µ2 (det F)2 F−1 : F−1 F−T +µ1 F−µ2 (det F)2 F−T F−1 F−T .
In this paper we consider some problems whose solutions can be computed in two
dimensions. In these cases and without loss of generality, B0 = B02D ×[0, 1], for a Lipschitz
domain B02D ⊂ R2 , and we denote X = (X2D , X3 ). In particular, we shall consider
plane-strain problems, in which ϕ(X) = (ϕ2D (X2D ), X3 ), for some planar deformation
mapping ϕ(X2D ). In the following, we will generally omit the superscript 2D to indicate
the two-dimensional quantities.
3. The HDG Formulation
We next define the approximations to the solutions of our non-linear elasticity problem
provided by the HDG method. We begin by introducing the notation associated to the
HDG methods and the finite dimensional spaces they use. Then we recall the notion
of the discontinuous Galerkin derivative [12] and, inspired in the minimization problem
defining the solutions of the non-linear elasticity problem, introduce a discrete version
of the minimization problem. Finally, we obtain the equations satisfied by the critical
points and show that they define the HDG method proposed by S.-C. Soon [1].
5
3.1. Notation
Here, we give the basic definitions needed for the formulation of the HDG methods.
Let Th be a conforming finite element mesh on B0 . Each element K ∈ Th is assumed
to be an open polygon (in two dimensions) or polyhedron (in three dimensions) with an
orientable boundary ∂K. Set ∂Th := {∂K : K ∈ Th }. We denote by n : ∂Th → Rd , the
function which restricted to ∂K is nothing but the outward unit normal. Let F denote
an arbitrary element face, let F(K) be the set of faces F of the element K ∈ Th , set
Γh := {F ∈ F(K) : K ∈ Th } and define
ΓIh = {F ∈ Γh : F ⊂ B0 },
Γdh = {F ∈ Γh : F ⊂ ∂d B0 },
ΓN
h = {F ∈ Γh : F ⊂ ∂N B0 }.
I
d
Note that Γh = ΓIh ∪ Γdh ∪ ΓN
h . Here Γh is the set of interior faces, Γh is the set of faces
N
lying on the Dirichlet boundary and Γh the set of faces lying on the Neumann boundary.
We use the following notation for the integrals:
Z
X
(a, b)K :=
ab dV
and
(a, b)Th :=
(a, b)K ,
K
K∈Th
Z
ha, bi∂K :=
ab dS
and
ha, bi∂Th :=
∂K
X
ha, bi∂K ,
K∈Th
Z
ha, biΓ :=
ab dS
for
Γ ⊆ ∂B0 .
Γ
Similar notation is used when Rthe integrands are vector-valued
or matrix-valued funcR
tions. For example, (a, b)K = K a · b dV and (A, B)K = K A : B dV .
b h ), in the
The HDG method seeks approximations to (∇ϕ|B0 , ϕ|B0 , ϕ|Γh ), (Dh , ϕh , ϕ
b h ) defined by
space (Gh , Vh , V
Gh := {G ∈ [L2 (B0 )]d×d : G|K ∈ [P k (K)]d×d ∀K ∈ Th }
Vh := {v ∈ [L2 (B0 )]d :
b h = {b
V
v ∈ [L2 (Γh )]d :
v|K ∈ [P k (K)]d ∀K ∈ Th }
b |F ∈ [P k (F )]d ∀F ∈ Γh }.
v
Here, P k (A) is the space of a polynomials of degree k defined on the set A. Note here
that functions in Gh and Vh are allowed to be discontinuous across element faces in ΓIh .
We can incorporate the Dirichlet boundary conditions into the space by introducing the
b h (ϕd ) := {b
b h : hb
b h }.
affine space V
v∈V
v, µiΓdh = hϕd , µiΓdh ∀µ ∈ V
3.2. The HDG Method as a Minimization Problem
To define the HDG method, we are going to use a discrete version of the minimization
problem defining the solutions of non-linear elasticity. To do that, we first have to note
that we are approximating the continuous displacement ϕ ∈ V by the discontinuous
b h . This has two main consequences. The first is that an
b h ) ∈ Vh × V
function (ϕh , ϕ
6
extra term needs to be added to the original functional in order to control the size of
the discontinuities of the approximation. The second, that we must come up with an
b h ). We define such
approximation of the gradient ∇ϕ defined solely in terms of (ϕh , ϕ
h bh
approximation as the element DDG (ϕ , ϕ ) of the space Gh determined as the solution
of
(6)
b h ), G)Th = −(ϕh , ∇ · G)Th + hb
ϕh , Gni∂Th
∀ G ∈ Gh .
(DDG (ϕh , ϕ
Note that DDG (ϕ, ϕ) is actually a good approximation of ∇ϕ. Indeed, DDG (ϕ, ϕ) is
nothing but the L2 -projection of ∇ϕ into the space Gh , as can be easily seen after a
simple integration by parts. This classical DG discretization of the gradient [35] has been
called the DG-derivative and has been used to define DG approximations for non-linear
elasticity in [12].
Next, let us introduce the discrete version of the potential energy functional I[ϕ]. To
emphasize its relation with its discrete counterpart, we begin by rewriting the functional
using our new notation as follows:
I[ϕ] := (W (∇ϕ), 1)B0 − (f , ϕ)B0 − hT, ϕi∂N B0
= (W (∇ϕ), 1)Th − (f , ϕ)Th − hT, ϕi∂N B0 .
ˆ h (g) → R is then defined by
The discrete energy functional Ih : Vh × V
b h ] :=(W (DDG (ϕh , ϕ
b h )), 1)Th − (f , ϕh )Th − hT, ϕ
b h i∂N B0
Ih [ϕh , ϕ
1
b h ), (ϕh − ϕ
b h )i∂Th .
+ h τ (ϕh − ϕ
2
(7)
b h ] as the sum over all the elements
Note that we can rewrite the functional Ih [ϕh , ϕ
b h ] defined by
K of the triangulation Th of similar, “elementwise” functionals IK [ϕh , ϕ
b h ] :=(W (DDG (ϕh , ϕ
b h )), 1)K − (f , ϕh )K − hT, ϕ
b h i∂N B0 ∩∂K
IK [ϕh , ϕ
1
b h ), (ϕh − ϕ
b h )i∂K .
+ h τ (ϕh − ϕ
2
Note also that the last term in the definition of the functional is meant to represent
an “energy” associated to the jumps of the approximation to the displacement across
interelement boundaries. For this reason, we choose the “stabilization” function τ to be
a scalar-valued strictly positive function which, for simplicity, is assumed to be constant
on each face of Γh .
b h ) to
We are now ready to define the HDG method. The approximation (Dh , ϕh , ϕ
b
the exact solution (∇ϕ|B0 , ϕ|B0 , ϕ|Γh ) is the element of the affine space (Gh , Vh , Vh (g))
b h ) and, (b) (ϕh , ϕ
b h ) is a local minimum of the
that satisfies that: (a) Dh = DDG (ϕh , ϕ
h
b h (g).
b ] on Vh × V
discrete potential energy functional Ih [ϕh , ϕ
3.3. The Euler-Lagrange Equations
Let us now find the equations satisfied by the critical points of Ih . The equations are
b h ], namely,
the Euler-Lagrange equations for the functional Ih [ϕh , ϕ
∂
b h (0),
b h, ]
b h ) ∈ Vh × V
Ih [ϕh, , ϕ
=0
∀ (δϕh , δ ϕ
∂
=0
7
b h, ) := (ϕh , ϕ
b h )+ (δϕh , δ ϕ
b h ). We solve these equations to find (ϕh , ϕ
b h ),
where (δϕh, , δ ϕ
through a Newton’s method described in §A. The Euler-Lagrange equations are contained in the following result.
b h (g)
b h ) of the functional Ih [ϕh , ϕ
b h ] over Vh × V
Theorem 1. The critical points (ϕh , ϕ
must satisfy that
∂W
b h )), DDG (v, µ)
b h , v − µi∂Th = (f , v)Th + hT, µi∂N B0
(DDG (ϕh , ϕ
+ τ h ϕh − ϕ
∂F
Th
(8)
b h (0). Alternatively, stated in a more standard HDG form,
for all (v, µ) ∈ Vh × V
h
b n, vi∂T = (f , v)T
(Ph , ∇v)Th − h P
h
h
h
b n, µi∂T = hT, µi∂ B
hP
N 0
h
∀ v ∈ Vh ,
(9a)
b h (0),
∀µ∈V
(9b)
b h n := Ph n − τ (ϕh − ϕ
b h ) on ∂Th , and Ph ∈ Gh satisfies that
where P
b h )), G)Th
(Ph , G)Th = (∂W/∂F(DDG (ϕh , ϕ
∀G ∈ Gh .
(9c)
b h ) is given by (6).
The DG-derivative DDG (ϕh , ϕ
Before proving this result, let us briefly comment on these equations. Equation (8)
states the balance of virtual work, as commonly found in numerical methods for nonlinear elasticity. The second set of equations coincides with the definition of the HDG
method proposed by S.-C. Soon [2, 1] for a slightly different definition of the stabilization function τ ; there τ is taken to be a fourth order tensor, whereas here we simply
take it to be a scalar. Equation (9a) can be interpreted as the static equilibrium over
every subset of elements in the mesh. Equation (9c) defines the approximation Ph to
the first Piola-Kirchhoff stress tensor. Finally, (9b) is nothing but a discrete version of
the transmission condition for the first Piola-Kirchhoff tensor, and, in the language of
discontinuous Galerkin methods for conservation laws, makes the method locally conserb h n single valued1 ; it also enforces the boundary
vative by rendering the numerical trace P
b h,
condition on ∂N B0 . To see this more clearly, note that, for any µ ∈ V
h
b n, µi∂T =
hP
h
X
K∈Th
h
b n, µi∂K =
hP
X
h,+ +
b
hP
b h,− n− , µiF + h P
b h n, µi∂ B ,
n +P
N 0
F ∈ΓIh
h,±
b
where P
are the traces from each of the sides of the face F and n± are the corresponding normals.
To prove Theorem 1 we need to obtain the first variation of the functional Ih , which
we do in the next result.
1 Only
b h n is defined, not that of P
b h.
the value of P
8
b h (g) and (δϕh , δ ϕ
b h ) be any element of Vh × V
b h ) any element of
Lemma 1. Let (ϕh , ϕ
b h (0). Then, we have that
Vh × V
∂W
∂
h, h
h
h,
h
h
b ]
b )), DDG (δϕ , δ ϕ
b )
=
Ih [ϕ , ϕ
(DDG (ϕ , ϕ
∂
∂F
=0
Th
b h , δϕh − δ ϕ
b h i∂Th − (f , δϕh )Th − hT, δ ϕ
b h i∂N B0
+ τ h ϕh − ϕ
b h n, δϕh i∂T −(f , δϕh )T +hP
b h n − T, δ ϕ
b h i∂Th ,
= (Ph , ∇δϕh )Th −h P
h
h
b h, ) := (ϕh , ϕ
b h ) + (δϕh , δ ϕ
b h ).
where (ϕh, , ϕ
Theorem 1 follows from this lemma by simply equating the first variation of Ih to zero
b h (0). The proof of Lemma 1 is included in
b h ) := (v, µ) ∈ Vh × V
and by taking (δϕh , δ ϕ
§B.
Remark. There is a close resemblance between the method for nonlinear elasticity discussed here and that in [28]. However, they are not identical. The method in [28] assumes
a customary additive decomposition of the strain energy density into an isochoric and a
volumetric part, namely, there it is assumed that
W (F) = U (det F) + Wiso (F)
for suitable functions U and Wiso , from where the pressure and the first Piola-Kirchhoff
stress tensor follow as
P(F, p) =
∂Wiso
(F) + p det F F−T ,
∂F
p = p(J) = U 0 (J).
(10)
The method in [28] then introduces an approximation of the pressure ph ∈ Vh = {v ∈
L2 (B0 ) : v|K ∈ P k (K) ∀K ∈ Th } found as
(ph , q)Th = (U 0 (det Fh ), q)
∀q ∈ Vh .
(11)
b h ). The equations of the
where to simplify the formula we identified Fh = DDG (δϕh , δ ϕ
h
method are then precisely (9a) and (9b) with P defined as
(Ph , G)Th =
∂Wiso
(Fh ) + ph det Fh F−T
,
G
h
∂F
Th
∀G ∈ Gh .
In contrast, for the method in this manuscript introducing (10) in (9c) leads to
(Ph , G)Th =
∂Wiso
(Fh ) + U 0 (det Fh ) det Fh F−T
,
G
h
∂F
Th
∀G ∈ Gh .
Clearly the definition of the approximate pressure through (11) does not guarantee that
6∈ P 0 (K) for each
the last two expressions are identical, essentially because det Fh F−T
h
element K.
9
0
20
40
60
80
100
120
140
160
180
200
0
50
100
150
200
(a) Coarsest mesh in some of the numerical (b) Sparsity structure of the global stiffness
matrix for the mesh in Fig. 1a.
h examples.
Figure 1: The coarsest mesh adopted for some of the numerical example is shown on the left. Other
meshes, especially for the determination of convergence rates, were constructed by sequentially and
uniformly subdividing this mesh. The sparsity structure for the stiffness matrix of the HDG method
with k = 1 on the mesh on the left is shown on the right. We also indicate the degree of freedoms
b h.
corresponding to ϕh and ϕ
4. A Non-Convergent Example and Its Relation to Stabilization
A weakness of non-conforming approximations such as DG or HDG for non-linear
elasticity problems is the possible loss of numerical stability, i.e., the appearance of
buckling modes that the exact solution does not display, if the stabilization parameter τ
is not large enough. Such loss of stability could lead to failure to converge as the mesh
is refined, as the example here shows.
Consider the simplest scenario in which the minimizer of the exact problem ϕ ∈ Vq
is unique, for Vq = W 1,q (B0 , R3 ) ∩ V and 1 ≤ q < ∞, and is such that there exist T > 0,
1 ≤ p ≤ q, and a continuous, monotonically increasing function f : [0, ∞) → [0, ∞) with
f (0) = 0 satisfying that
kΨ − ϕkVp < f (I[Ψ] − I[ϕ])
(12)
for all Ψ ∈ UT = I −1 [ [I[ϕ], I[ϕ] + T ) ] ∩ Vq . Finally, we assume that ϕd is affine, so
that Vh ∩ Vq is non-empty.
Consider then a standard conforming finite element method, which seeks
ϕhc = arg
max
Ψ∈(Vq ∩Vh )
I[Ψ].
(13)
Additionally, assume that the method is such that for a sequence of meshes {Th }h with
˜ h }h such that
h → 0+ , there exists a sequence {Ψ
˜ h ] = I[ϕ].
lim I[Ψ
h→0+
(14)
This condition guarantees that the finite element spaces have the ability to approximate
the minimum value of the functional. In this case, the method inherits the trivial stability
10
estimate
kΨh − ϕkVp < f (I[Ψh ] − I[ϕ])
for all Ψh ∈ UT ∩ Vh . Because of (14), for any 0 < t ≤ T the set Ut ∩ Vh is non-empty
for h small enough, and hence ϕhc ∈ Ut . It then follows that kϕhc − ϕkVp < f (t), a
convergence result.
In the non-conforming, HDG method, such stability estimate is not necessarily inherited. To see this, first notice that
I[Ψ] = Ih [Ψ, Ψ|∂Th ]
b h ) such that
for any Ψ ∈ Vq . The problem is, essentially, that there may be (Ψh , Ψ
b h ) − (ϕ, ϕ|∂T )k
k(Ψh , Ψ
>C>0
b
h
Vh ×V(0)
(15)
b h (0), and
in a suitable defined norm for Vh × V
b h ] ≤ Ih [ϕ, ϕ|∂T ] = I[ϕ]
Ih [Ψh , Ψ
h
(16)
for all h small enough. In other words, the presence of discontinuities may enable Ih to
attain values that are even smaller than that of the exact solution, on deformations that
are far from it regardless of the value of h.
We show next, with a numerical example, that this can actually happen. Consider a
body with B0 = [0, 1]2 and a strain-energy density function of the form
c (C) =
W
λ
µ1
µ2
C: I +
(C − I) : (C − I)
+
1/2
2
8
(det C)
(17)
where I is the 3 × 3 identity matrix, and µ1 , µ2 > 0. The material constants are λ = 0.31
and µ1 = 0.01, µ2 = 0.05. The body is fixed all along ∂B0 , namely, ϕ(X) = X on ∂B0 ,
and we set f = 0 and T = 0.
Lemma 2. For the problem described above, the unique minimizer of (2) in V4 = V ∩
W 1,4 (B0 , R3 ) √
is the deformation ϕ(X) = X. Furthermore, if λ > µ1 , it satisfies (12)
with f (t) = C t and p = 2, for some constant C > 0.
Before proceeding to the proof, notice that, as argued above, for this problem a conform˜ h (X) = X ∈ Vh would be stable, since (14) would
ing method given by (13) for which Ψ
be trivially satisfied.
proof. Recall first that for two-dimensional problems we have assumed plane strain
conditions, and such deformations have the form Ψ = (Ψ2D (X1 , X2 ), X3 ). Therefore, in
the Cartesian basis parallel to the coordinate axes, the tensor C has the form


C11 C12 0
[C] = C12 C22 0 ,
0
0
1
p
√
2 >
where C11 > 0, C22 > 0, and for admissible deformations, J = det C = C11 C22 − C12
0 a.e. in B0 .
11
The proof proceeds by bounding the difference in potential energies between ϕ and
any Ψ ∈ V4 from below, to obtain (12), from where the uniqueness of the minimizer
follows directly. Let C = ∇ΨT ∇Ψ, and recall that ∇ϕ = I in B0 . Then,
µ2
I[Ψ] − I[ϕ] ≥
8
1
− 1 dV
(C − I) : (C − I) dV + (λ − µ1 )
(det C)1/2
B0
B0
Z 1
−1/2
+ µ1
(C − I) : I + (det C)
− 1 dV
2
B0
Z
Z
(18)
The second and third terms on the left hand side are non-negative, as we describe next.
The non-negativity of the second term follows from noticing that [see e.g. 36, §4.3.2, Thm
3.2]
Z
Z
(det C)1/2 dV =
B0
det ∇Ψ dV = |B0 |,
B0
and hence, from Jensen’s inequality applied to the convex function z(x) = 1/x for x > 0,
that
Z
1
1
1
= 1.
dV ≥ 1 R
|B0 | B0 (det C)1/2
(det
C)1/2 dV
|B0 | B0
To see that the third term is non-negative, notice that
Z
Z
1
1
−1/2
(C − I) : I + (det C)
≥
(C11 + C22 − 2) + (C11 C22 )−1/2 ,
2
2
B0
B0
and that
(C11 + C22 − 2) + (C11 C22 )−1/2 ≥ 1
for C11 ≥ 0 and C22 ≥ 0. We can then conclude that
I[Ψ] − I[ϕ] ≥
µ2
kC − Ik2L2 (B0 ) ≥ C1 kF − Ik2L2 (B0 ) ≥ C2 kΨ − ϕk2H 1 (B0 ) ,
8
(19)
for all Ψ ∈ V4 . The positive constants C1 and C2 follow from a non-linear Korn inequality
shown in [37, Thm. 1], and Poincar´e inequality in W 1,2 (B0 , R3 ), respectively.
Next, we look at what happens when the solution of this problem is approximated
with the HDG method. The exact solution of this problem belongs to the HDG space,
b h (X). However, when τ is not large enough, it is not the
namely, (X, X|∂Th ) ∈ Vh × V
minimizer of Ih . It is only a stationary point.
Lemma 3. The exact solution (X, X|∂Th ) satisfies the Euler-Lagrange equations (8),
and is thus a stationary point of Ih .
proof. In this case DDG (X, X|∂Th ) = I in Th , which follows from (25), and ∂W/∂F(I) =
(µ1 − λ)I. Thus, since f = 0 and T = 0, (8) reads
((µ1 − λ)I, DDG (v, µ))Th = 0
b h (0). This equation is trivially satisfied, as a consequence of the
for all (v, µ) ∈ Vh × V
definition of the DG-derivative in (6).
12
Figure 2: Equilibrium configurations for Ih and τ = 1, with a lower potential energy than the exact
solution, for two different meshes. Contours of the modulus of the displacement field are shown.
We show next that (X, X|∂Th ) is not a global minimizer. We do so by both computing
other stationary points of Ih with smaller energy, and by plotting the value of Ih nearby
the exact solution. From a mechanical perpective, the exact solution is thus an unstable
equilibrium, and small perturbations induce large displacements, or numerical buckling.
We computed values of Ih and corresponding minimizers by letting k = 1 and τ = 1,
on meshes obtained as uniform, self-similar refinements of the mesh in Fig. 1, starting
from the first refinement. By starting Newton’s iterations from an initial guess distinct
from the exact solution, we computed numerical solutions that differ from it. Two of
such solutions for two different mesh sizes are shown in Fig. 2. It is evident from the lack
of four-fold symmetry, that other equilibrium configurations exists, just as with buckling
modes.
b h ) to the exact solution was measured in the norm
The distance or error of (ϕh , ϕ
b h )kVh ×V(0)
b h )kL2 (B0 )
k(ϕh , ϕ
= kϕh kL2 (B0 ) + kDDG (ϕh , ϕ
b
as a function of the mesh size h. The results are shown in Fig. 3a. In the range of mesh
sizes explored, there is no indication of the error decreasing as h decreases. In fact, the
figure suggests that the error may remain away from zero regardless of the value of h. In
other words, in addition to the exact solution, there exist at least one more stationary
point of Ih . Fig. 3b shows the energy of each one of the numerical solutions in Fig. 3a.
As a reference, the value of the potential energy at the exact solution is also shown. The
figure clearly demonstrates that for each mesh there exists at least another stationary
point for which Ih takes smaller values than at the exact solution. Together, both figures
suggest that the scenario hypothesized in (15) and (16) is realized in this example.
bh)
To explore the values of the energy around the exact solution for each h, let (uh , u
be the eigenmode corresponding to the minimum eigenvalue of the stiffness matrix at the
exact solution (X, X|∂Th ), normalized so that the sum of the squares of all nodal values of
b h ) is 1. For all the meshes in this study, these eigenvalues are consistently negative,
(uh , u
decreasing from ≈ −0.6 to ≈ −0.9 across five subdivisions of the initial mesh. Fig. 4a
13
0.475
0.32
0.47
0.315
0.465
0.31
HDG solution
I [X] = I h [(X, X|∂T h )]
Ih
b h ) − I|| L 2 (B 0 )
||ϕ h − X|| L 2 (B 0 ) + ||D DG (ϕ h , ϕ
b )] as a function of α, for different values of h.
shows the value of Ih [(ϕ, ϕ|∂Th ) + α(u, u
As a reference, the value of the potential energy at the exact solution is also shown. As
expected from Lemma 3, the functional Ih is essentially constant near α = 0. However,
for larger values of α, the value of Ih is smaller than its value at the exact solution,
regardless of the value of h. The use of a large enough stabilization parameter eliminates
these “spurious” minimizers, as Fig. 4b shows.
0.305
0.46
0.455
0.3
0.45
0.295
0.445 −3
10
−2
−1
10
0.29 −3
10
0
10
10
−2
−1
10
0
10
h
10
h
(a) Error in the numerical solution with re- (b) Value of Ih at each one of the deformations in
spect to the exact solution, the identity map, Fig. 3a. The value of Ih at the exact solution is
as a function of the mesh size.
shown as well.
Figure 3: The two figures together suggest the realization of (15) and (16), namely, the existence of
minimizers of Ih that remain away from the exact solution regardless of the value of h. Thus, the exact
solution becomes a saddle point for Ih , instead of a local or global minimizer, and it is hence unstable.
Setting the value of τ large enough removes this instability, as shown in Fig. 4b.
0.33
1
0.9
0.32
0.8
τ = 1e − 02
τ =1
τ = 10
τ = 100
I [X] = I [(X, X|∂T h )]
0.31
Ih
Ih
0.7
0.3
0.6
0.5
0.29
0.28
0.27 −3
10
h = 0.2042
h = 0.1021
h = 0.051031
h = 0.0255
I [X] = I h [(X, X|∂T h )]
0.4
0.3
−2
−1
10
10
α
(a) Case τ = 1.
0
10
0.2 −3
10
−2
−1
10
10
0
10
α
(b) Case h ≈ 0.05
b h )] as a function of α. Here (uh , u
b h ) is a disFigure 4: Value of the energy Ih [(X, X|∂Th ) + α(uh , u
placement field in the direction in which the energy Ih decreases the most, namely, the eigenmode
corresponding to the minimum eigenvalue of the stiffness matrix for each h. The existence of a minimum
away from α = 0, or the exact solution, is clearly shown in Fig. 4a. Such numerical instabilities are
removed by increasing the value of τ , as shown in Fig. 4b.
To summarize, this example shows that the HDG method does not necessarily inherit
the stability properties of the problem being solved, and that a large enough value of
τ needs to be found for this to be the case. In some situations, such as not too close
14
to buckling loads, the adaptive stabilization strategy in [21, 22] could likely be used
to automatically adjust the value of the stabilization parameter as the load progresses.
We do not apply such strategy here. Instead, the values of the stabilization parameters
throughout this paper have been the result of a few trials.
Let us point out that the stabilization function on a given face used in [21, 22] was
τ = β/h, where h is the smallest diameter of the elements sharing the face in question and
β is a parameter scaling like the elastic moduli. If we adopt this choice of stabilization
function, the optimal convergence of the stresses is lost. To avoid this unwanted behavior,
τ must be independent of h; see the numerical results in [2, 1] and ours in §5. In other
words, we should take τ = β/L, where β scales like the elastic moduli and L like a length.
At this point, however, how to pick the best values of β and L remains an open problem.
5. Numerical Examples
Next, we showcase the performance of the proposed method for various elasticity problems. We first examine examples in two dimensions, with two different non-quadratic
strain-energy density functions leading to non-linear behavior. In these examples we
investigate the numerically obtained convergence rates of the method, as well as its
behavior under locking situations due to kinematic constraints imposed by near incompressible behavior. We additionally showcase the deformations obtained under very large
strains, including an example in two-dimensional cavitation. Two three-dimensional examples demonstrate deformations in which controlled buckling of the structure occurs.
In a simpler example, we then compare the performance of the proposed HDG method
against a standard discontinuous Galerkin and a standard continuous Galerkin method,
and obtain rather encouraging results for the types of low order elements we adopted.
All of the forthcoming numerical examples were performed by setting k = 1, namely,
b h . Also, all stress contours were
elements with affine shape functions for both ϕh and ϕ
h bh
plotted by interpolating P(DDG (ϕ , ϕ )) at the mesh nodes, instead of by computing
Ph . Because the difference between the two is quite small, no explicit distinction was
made in the figures.
As expected, the approximations for the two-dimensional case display similar properties than those obtained in [1], although the example problems are rather different.
The three-dimensional numerical experiments are obtained here for the first time for this
method.
5.1. Convergence Rates
We first illustrate the numerically observed convergence rates for a simple plane-strain
problem in non-linear elasticity for which we construct an analytical solution, and for
two different non-linear materials. To this end we consider families of structured meshes
over B0 = (1, 2) × (1, 2), as depicted in Figure 1a.
To construct problems with an analytical solution, we first choose the deformation
mapping, and then compute the corresponding body forces needed for it to be a solution.
Additionally, we set ∂d B0 = ∂B and evaluate ϕd on ∂d B0 to coincide with the values of
the exact solution. The chosen exact solution is ϕ(X, Y ) = (X 2 , Y 2 ), where (X, Y ) are
the Cartesian coordinates in R2 . For it to be a solution, the necessary body forces for a
15
Mooney-Rivlin material are given by
f (X, Y ) =
λ [1 − log(4XY )] + µ1 + 4X 2 (µ1 + µ2 )
eX
2X 2
λ [1 − log(4XY )] + µ1 + 4Y 2 (µ1 + µ2 )
+
eY , (20)
2Y 2
where eX and eY are the unit vectors in the X and Y directions, respectively. The
corresponding body forces for a neo-Hookean model follow by setting µ2 = 0. In the
following computations we set λ = 1 and µ1 = 1 for both material models.
We tested the proposed HDG formulation with τ = 1 and k = 1, namely, affine
b h . By sequentially subdividing each triangle of the mesh
polynomials for both ϕh and ϕ
into four similar ones, beginning with the mesh in Fig. 1a, we obtained a family of
meshes with uniformly decreasing mesh size h. The errors ||ϕ − ϕh ||L2 (B0 ) and ||∇ϕ −
b h )||L2 (B0 ) as a function of the mesh size h are shown in Fig. 5, while the
DDG (ϕh , ϕ
numerical values and the computed convergence rates are shown in Tables 1 and 2. For
b h ) converge with order k + 1 = 2. The
both material models, both ϕh and DDG (ϕh , ϕ
convergence of the stress tensor P is that of DDG . This follows from the fact that P
is a uniformly continuous function over the closure of an open neighborhood of the set
{F ∈ R2×2 | F = ∇ϕ(X, Y ), (X, Y ) ∈ B0 }.
0
−1
10
10
−1
−2
10
10
−2
−3
10
10
−3
−4
10
10
−4
−5
10
10
−5
10 −2
10
−6
−1
10
(a) Neo-Hookean material
10
−2
10
−1
10
(b) Mooney-Rivlin material
Figure 5: Convergence curves towards an analytical, manufactured solution, for both material models.
The expected convergence rates for an HDG method are observed.
5.2. Incompressible Limit for Large Deformations
We next explore the performance of the method under near incompressibility conditions. HDG methods do not display locking in near incompressible situations for problems
in isotropic linear elasticity [1]. As we show next, locking does not appear either for the
non-linear elastic problem we consider here, which involves quite large deformations of
the body.
The fundamental difference between the linear and non-linear elasticity cases is that
the incompressibility constraint in the latter is non-linear. More precisely, in a perfectly
incompressible situation we should have det F = 1 almost everywhere in the domain,
instead of a divergence-free condition.
16
Table 1: Convergence rates for a neo-Hookean material.
h0 /h
1
2
4
8
16
||ϕ − ϕh ||L2 (B0 )
Error
Order
2.380e-02
5.893e-03
1.467e-03
3.661e-04
9.142e-05
2.013
2.005
2.003
2.001
b h )||L2 (B0 )
||∇ϕ − DDG (ϕh , ϕ
Error
Order
1.673e-02
4.625e-02
1.256e-02
3.370e-03
8.964e-04
1.855
1.880
1.897
1.910
Table 2: Convergence rates for a Mooney-Rivlin material.
h0 /h
1
2
4
8
16
32
||ϕ − ϕh ||L2 (B0 )
Error
Order
6.598e-03
1.648e-03
4.092e-04
1.021e-04
2.549e-05
6.368e-06
2.005
2.005
2.003
2.017
2.000
b h )||L2 (B0 )
||∇ϕ − DDG (ϕh , ϕ
Error
Order
3.598e-02
9.468e-03
2.491e-03
6.539e-04
1.698e-04
4.300e-05
1.926
1.925
1.929
1.944
1.981
The test problem consists of an annulus made of a nearly incompressible neo-Hookean
material. Let (R, Θ) denote the polar coordinates in the reference configuration. Then,
the reference configuration B0 is the set of points for which 0 < R0 = 1/2 < R < R1 = 1.
The annulus is deformed by setting ϕd (X) = Xr0 /R0 on ∂d B0 = SR0 , and T = 0 on
∂N B0 = SR1 , where SR is the circumference of radius R centered at the origin.
In the perfectly incompressible case, this problem has a simple analytical solution, see
e.g. [12]. In the near incompressible case an analytical solution seems to be unavailable.
However, assuming that the solution is axisymmetric, it is possible to reduce the solution
of the problem to solving an ordinary differential equation. Its solution can be computed
very accurately with a number of standard algorithms. This is the approach we adopted
here to be able to construct what we term a “quasi-analytical” solution, which we subsequently use as the reference solution to evaluate the convergence rate of the proposed
HDG method.
We detail the reduction of the problem to an ordinary differential equation next.
First, we assume that the solution is axisymmetric, namely, that ϕ(X) = ϕ(|X|)X/|X|
for some function ϕ : [R0 , R1 ] → R+ . The deformation gradient for such deformation
mapping is
ϕ(|X|)
∇ϕ(X) = ϕ0 (|X|)eR ⊗ eR +
eΘ ⊗ eΘ ,
|X|
where (eR (|X|), eΘ (|X|)) is the standard orthonormal basis at point X in polar coordinates. Under these deformations, the first Piola-Kirchhoff stress tensor for a neo-Hookean
material takes the form
P = PRR eR ⊗ eR + PΘΘ eΘ ⊗ eΘ ,
17
(21)
Figure 6: Reference (center) and deformed configurations, in scale, for a mesh with 14628 triangles and
λ = 1666.44. The inset shows an enlargement of a region of the deformed annulus. Contour levels of
h are shown.
PRR
where
λ log
PRR
=
Pθθ
=
ϕ(R)ϕ0 (R)
R
1
0
+
ϕ
(R)
ϕ0 (R)
ϕ0 (R)
0
(R)
λ log ϕ(R)ϕ
R
ϕ(R)
R
+ µ1 −
+
.
ϕ(R)/R
ϕ(R)
R
+ µ1 −
The ordinary differential equation on ϕ follows from imposing the equilibrium condition
(or one of the Euler-Lagrange equations for (3)) ∇ · P = 0 :
dPRR
PRR − Pθθ
+
= 0.
dR
R
(22)
This equation needs to be solved with the boundary conditions ϕ(R0 ) = r0 and PRR (R1 ) =
0 for each value of the elastic moduli that we test. In all cases we obtained quasi-analytical
solutions with an accuracy of at least 10 significant digits for ϕ and ϕ0 .
We next proceed to compare the results from the quasi-analytical solution with those
obtained from the HDG method, as λ → ∞, which is the near incompressibility limit
for this material model. In the examples below we set r0 = 1.5, τ = 1, µ1 = 0.333 and
k = 1. Because we are imposing the deformation mapping on ∂d B0 , locking would be
manifested as an unbounded increase of the stress field as λ → ∞ for a fixed mesh.
Results of this numerical experiments are shown in Figs. 6, 7, and 8. The reference
and deformed configuration are shown in scale in Fig. 6, together with a detailed view
h
of a part of the ring. Shown as well are the contours of PRR
, which are quite smooth,
and as we show later, converge to the exact solution regardless of the value of λ; there
18
b h ) over the reference configuration, for meshes with 1182
Figure 7: Contour plots of det DDG (ϕh , ϕ
(left) and 3792 (right) elements and λ = 1666.44. It should be close to 1 everywhere, and it ranges from
0.998 to 1.
−1
10
||ϕ − ϕh || L 2 (B 0 )
||ϕ − ϕh || L 2 (B 0 )
b h )|| L 2 (B 0 )
||∇ϕ − D D G(ϕh , ϕ
b h )|| L 2 (B 0 )
||∇ϕ − D DG (ϕh , ϕ
−2
10
||P − P h || L 2 (B 0 )
||P − P h || L 2 (B 0 )
−1
10
−3
−4
10
Error
Error
10
O(h 2 )
−2
10
−5
10
−6
10
−7
10
−3
−1
10
10
2
10
h
3
10
4
5
10
10
6
10
7
10
λ
(a) Error as a function of h, for λ = 1666.44.
(b) Error as a function of λ, for h = 0.1.
Figure 8: The HDG method converges with the expected rates in a near incompressible situation (left).
Such convergence is apparently uniform on λ, as the figure on the right shows. This is consistent with
a locking-free behavior of the method even in this highly non-linear case. The error is computed as the
difference between the HDG solution and the quasi-analytical solution.
19
b h on two different meshes over
is no sign of locking. Contour plots for det DDG (ϕh , ϕ
the reference configuration are shown in Fig. 7, for λ = 1666.44.
Convergence curves are shown in Fig. 8, as well as evidence consistent with such
convergence being uniform with respect to λ, or in other words, locking-free. The error
is computed as the difference between the HDG solution and the quasi-analytical solution.
As mentioned earlier, in this example locking would be manifested as a steady growth of
Ph with λ for the same mesh. This figure shows that this is not the case.
The solution of problems for different values of r0 was possible up to values of
r0 = 2.48. For larger values, we found difficulties for the Newton-Raphson iterations
to converge, so an alternative root-finding strategy might be needed then.
5.3. Examples with Large Deformations in Two- and Three-dimensional Geometries
In the following sections we showcase examples in two- and three-dimensional geometries in which the elastic body undergoes very large deformations. The goal of these
examples is to heuristically illustrate that the method works well in rather demanding
situations, and that the numerical solutions could particularly benefit from the discontinuities allowed in the deformation mapping.
5.3.1. 2D: Twist of an Annulus
The first example explores the performance of the method under large shear strains.
To this end, we consider a circular annulus of outer radius 1 and inner radius 0.333 made
of a neo-Hookean material with λ = 1 and µ = 1, see Fig. 9. The inner boundary of
the annulus is kept fixed, and the outer one is quasi-statically rotated around the center
of the annulus in the counterclockwise direction. We apply this load through Dirichlet
boundary conditions on both boundaries, rotating the outer boundary in angle steps of
π/60. The stabilization parameter τ was taken to be equal to 10.
This example is difficult to perform because of the large shear strains that appear
near the inner boundary, see Fig. 9. There, large discontinuities appear across elements,
see Fig. 9c. To verify the accuracy of the computations, we performed simulations for
a mesh, and for one uniform subdivision of such mesh, and computed a quasi-analytical
solution, similarly to the approach in §5.2, exploiting the axi-symmetry of the problem.
We skip the details of the computation here. The values of PRθ and PθR computed from
the quasi-analytical solution and from the numerical solutions at all nodes of the two
meshes are shown in Fig. 10. A trend of convergence is observed when the mesh is
refined, with larger errors near the inner boundary, as expected. This suggests that the
numerical solution is meaningful, despite the large discontinuities near the boundary.
The presence of such discontinuities enables the simulation to continue for much
larger angles than with CG. The CG simulation for these meshes could not go pass a
105◦ rotation, since elements near the inner boundary would invert in such case. Instead,
issues with HDG appeared at approximately 185◦ rotation.
5.3.2. Three-dimensional Examples
Thus far all simulations presented have been restricted to two-dimensional examples,
for their lower computational cost makes them more convenient to study the numerical
properties of the method. We now illustrate the behavior of the method through a couple
of three-dimensional examples.
20
(a) After rotation by
π
.
2
(c) Closer look at ϕh
(b) After rotation by π.
Figure 9: Two snapshots of the deformed configuration. The color contours correspond to the norm
of the deformation, kϕh k. There are discontinuities throughout the domain, but those near the inner
boundary are particularly larger. Figure 9c shows the enlargement of the region marked in 9b for ϕh .
The discontinuity of the deformation mapping can be observed.
1.4
3
4728 elements
1182 elements
Quasi-analytical
1.2
4728 elements
1182 elements
Quasi-analytical
2.8
2.6
1
2.4
2.2
h
PθR
h
PRθ
0.8
0.6
2
1.8
0.4
1.6
0.2
1.4
0
−0.2
1.2
0.4
0.5
0.6
0.7
0.8
0.9
1
1
0.4
0.5
0.6
0.7
0.8
0.9
1
R
R
Figure 10: Comparison of two components of Ph among a quasi-analytical solution computed by exploiting the axi-symmetry of the problem, and the HDG solution on two different meshes. All nodal
values for the stress components are plotted. A convergence trend toward the quasi-analytical solution
as the mesh is refined can be inferred.
21
Figure 11: Snapshots of the deformed configuration of a hollow cylinder are shown, where the only
elements shown are those that are visible once all elements on the other half of the cylinder have been
removed. The inset s hows an enlargement of a region of the deformed cylinder, in which some of the
discontinuities in the deformation mapping can be observed. The contours represents ||ϕh − X||.
The first example consists of a hollow cylinder made of a neo-Hookean material deforming under compression and shear, with λ = 1 and µ = 0.1. The inner radius, outer
radius, and height of the cylinder are 0.75, 1, and 4, respectively. The mesh consists of
22,000 tetrahedral elements. All computations were performed with k = 1 and τ = 100.
The cylinder is deformed such that its bottom face is held fixed, while the top face is
translated both vertically and horizontally at a uniform rate per loading step, up to a
maximum displacement of 1.5 in both directions. The resulting deformations in the cylinder can be roughly described as compressive with shear, as revealed by the snapshots in
figure 11. These show the reference configuration as well as the deformed configuration
of the cylinder at three different instants of the loading path. Notice that the solution
displays an emerging localization of the deformation. The contours of a component of
the first Piola-Kirchhoff stress tensor at the same snapshots are shown in Fig. 12, all
drawn over the reference configuration.
In the second example we considered a hollow cylinder that has half the height of
that in the previous example, and is made of a neo-Hookean material with λ = 0.4,
22
-6
-18
-30
-42
Figure 12: Contour plots of n · Ph n over the reference configuration, for the same snapshots in Fig. 11.
Here n is the direction in which the cylinder is sheared.
Figure 13: Snapshots of deformed configurations of the hollow cylinder under torsional load. The last
snapshot shows the final deformation, in which the top face was rotated by 180◦ . The contours represent
||ϕh − X||.
µ = 0.4. We torsioned the cylinder by keeping the lower base fixed, and by imposing a
rigid rotation in the counterclockwise direction to the top base, for progressively larger
torsion angles; in both cases through Dirichlet boundary conditions. The maximum
rotation angle is 180◦ , and it is reached in 60 steps of 3◦ each. The rest of the boundary
remains traction free. The mesh is made of 22,000 tetrahedral elements, and we set k = 1
and τ = 100. Figure 13 shows snapshots of three deformed configurations of the cylinder.
Fig. 14 displays the stress contours over the reference configuration. In particular, we
h
show the value of the nodal interpolant of σθz
. This is the approximation to the shear
h
stress σθz , or the θz-component of the Cauchy stress tensor. The value of σθz
is computed
h
−T
h bh
by replacing F by DDG (ϕ , ϕ ) and P by P in σθz = det F eθ · PF ez , where ez
and eθ are the unit vectors in the axial and azimuthal directions, respectively. Notice
that this example includes large shear deformations, especially near the lower and upper
surface of the cylinder.
23
10
3.3
-3.3
-10
h , or the azimuthal-axial shear stress component of the Cauchy stress
Figure 14: Contour plots of σθz
tensor, over the reference configuration, for the snapshots in Fig. 13.
5.4. Cavitation
The last example involves the computation of a cavitating void. Cavitation refers
to a phenomenon observed in elastomers (see, e.g. [38]) in which, for sufficiently large
tensile stresses, small voids enlarge by one or more orders of magnitude. The numerical computation of cavitation is a difficult problem because the large expansion induces
very large strains on elements near the cavity. These difficulties were nicely illustrated
and discussed by Xu and Henao [39]. In particular, they highlight advantages of some
non-H 1 -conforming methods for cavitation problems, although H 1 -conforming isoparametric quadratic elements have also been found useful for these problems [40]. Threedimensional simulations with highly-refined meshes of 8-nodded bilinear hexahedra have
been reported as well [41].
In the following we showcase the performance of the method proposed here on two
cavitation examples in two-dimensions, and one in three-dimensions. For cavitation to
occur, a slower growth of the strain-energy density function with the deformation gradient
is needed. In this case, we adopt the model
W (F) =
3
2µ1
(F : F) 4 + U (det F),
5/4
3
(23)
where again µ1 is a material constant, and U is given in (5).
The two two-dimensional examples consist of a disc of unit radius with one or two
holes as the reference configurations. In the first example, a single hole of radius 0.2
is placed centered at (0.3, 0), while in the second example two holes of radii 0.22 and
0.2 are placed centered at (−0.3, 0) and (0.3, 0), respectively. In both cases, the discs
are deformed by leaving the walls of the holes traction-free, while imposing Dirichlet
boundary conditions ϕd (X) = αX on the outer wall |X| = 1, for progressively larger
values of α ≥ 1. In both cases, a mesh with 8987 elements with k = 1 was adopted, and
we set τ = 50, µ1 = 1, and λ = 1. Figures 15 and 16 show the results of the numerical
experiments for a value of α = 4.7 and α = 6.16, respectively. The deformations in
both cases are large, and the elements are very stretched. To reach these values of α, we
computed a total of 150 and 197 loading steps, respectively, by increasing the value of α
by 0.0313 on each step, and by adopting the solution of the last step as the initial guess
24
32
28
24
20
16
12
8
4
0
Figure 15: The reference configuration and the corresponding mesh are shown on the left, together with
h . Deformation of the first two-dimensional example for α = 6.16, on the right.
the contour plot for Pθθ
An enlargement of a region of the deformed configuration is shown at the center of the cavity.
for the solution the new step. Convergence of the non-linear solver further than these
values of α was rather delicate, and hence we did not proceed any further.
The reference configuration for the three-dimensional example is a sphere of radius
1 with two spherical cavities inside. The first cavity has radius 0.2 and is centered at
(−0.7, 0.7, 0) in a Cartesian coordinate system with origin at the center of the sphere. The
second cavity has radius 0.45, and is centered at (0.25, 0.25, 0.25). As in two-dimensions,
the sphere is deformed by leaving the walls of the holes traction-free, while imposing
Dirichlet boundary conditions ϕd (X) = αX on the outer wall |X| = 1, for progressively
larger values of α ≥ 1.
We performed simulations on a mesh with 47073 tetrahedra, see Fig. 17, and by
setting λ = 1, µ1 = 1, k = 1, and τ = 100. In this case, the value of α was increased by
0.0475, and difficulties to find the solution of the non-linear system appeared for values
of α greater than 2.85. Snapshots of the deformed configuration for α = 2.85 are shown
in Fig. 18. Contour plots of ||ϕ − X|| are shown as well.
We have yet to run convergence studies on these examples. However, the results shown
are very promising, and display the kind of added flexibility that non-H 1 conforming
methods, and in particular HDG, have for cavitation problems.
5.5. Efficiency
In this section we compare the relative efficiency of the HDG, DG and the CG methods
by plotting error versus the corresponding computational complexity. We take computational complexity to be proportional to the size of the stiffness matrix or the number of
non-zero entries of the stiffness matrix. Additionally, we also compare the error of these
methods as a function of the time for solving each non-linear system of equations. We
show results for two- and three-dimensional problems.
25
32
28
24
20
16
12
8
4
0
Figure 16: The reference configuration and the corresponding mesh are shown on the left, together with
h . Deformation of the second two-dimensional example for α = 4.70, on the right.
the contour plot for Pθθ
An enlargement of a region of the deformed configuration is shown at the center of the cavity.
Figure 17: Mesh of 47073 tetrahedra for the reference configuration of the sphere with two cavities.
26
Figure 18: Deformed configuration for the sphere with two cavities of Fig. 17 after stretching its outer
boundary to a radius nearly three times its original size (α = 2.85). The right figure shows the elements
that are visible once all elements on the other half of the sphere have been removed, while the left figure
shows the intersection of the sphere’s deformed configuration with two orthogonal planes through its
center. The original cavities in the interior of the sphere expanded even more than the sphere, and the
formation of a thin wall between the two is starting to appear, as in the two-dimensional case. The
contours for ||ϕ − X|| are shown as well.
We compared the three methods by solving the problem described in §5.1 for a
neo-Hookean material in two-dimensions. We formulated a similar problem in threedimensions, to obtain an exact solution of the form ϕ(X, Y, Z) = (X 2 , Y 2 , Z 2 ), where
(X, Y, Z) are Cartesian coordinates, and the reference configuration is the cube B0 =
(1, 2)3 . In two dimensions we adopted the meshes obtained from that in Fig. 1a and
its self-similar subdivisions. In three-dimensions, we chose an analog structured mesh of
tetrahedra and its subdivisions according to [42].
The size of the stiffness matrix for HDG is the sum of the number of degrees of
b h over the entire mesh, for k = 1. Notice that this choice is
freedom for ϕh and ϕ
the worst possible scenario for HDG, since it does not reflect the fact that the stiffness
matrix shows an essentially block-diagonal structure for the degrees of freedom of the
same elements. This structure enables the “static condensation” of the degrees of freedom
of ϕh , or alternatively, a very fast solution of these degrees of freedom. The CG method is
obtained in a standard way by adopting polynomials that are affine inside each element,
and are continuous across element boundaries. The size of the stiffness matrix in this
case is approximately equal to two or three times the number of nodes in the mesh,
depending on the spatial dimension. Finally, the DG method we use is the one proposed
in [12], with k = 1 as well; the size of its stiffness matrix is approximately 3 or 4 times
the number of elements in the mesh, depending on the dimension. The solutions with
DG were computed only for the two-dimensional cases.
The comparison are shown in Figs. 19, 20, 21 and 22. In two-dimensions, we computed with 6 meshes, obtained by uniform subdivisions of the initial mesh, for DG and
27
0
10
0
10
−1
10
−1
10
−2
−2
Error
Error
10
−3
10
−4
10
10
−3
10
−4
10
HDG
HDG
CG
10
CG
DG
−5
1
10
DG
−5
2
10
3
10
4
10
5
10
Size of the system
6
10
10
1
10
2
10
3
10
4
10
5
10
6
10
7
10
Non-zero entries in matrix
Figure 19: ||ϕ − ϕh ||L2 (B0 ) versus computational cost for a two-dimensional geometry.
HDG, and with 8 for CG. In three-dimensions, we computed with 4 meshes, again subdivisions of the original mesh. As seen in Figs. 19 and 21, the accuracy of the deformation
mapping ϕh for the same computational cost is better for CG throughout. The accuracy
of the approximation to ∇ϕ is much better in the case of HDG. In fact, in two-dimensions
the HDG solution is better just after the first subdivision, and in three-dimensions the
HDG solution is more efficient after around 40,000 degrees of freedom.
We compared the processor time consumed by the program to solve the non-linear
equations using Newton’s method with CG, HDG, and DG methods in two-dimensions,
and with the first two in three-dimensions. Comparisons of time are always fraught with
issues related to how good the implementations and optimizations of each code are. Our
perspective here is to see whether the computing times were at least in the same order
of magnitude, rather than more accurate performance information. We believe that the
results shown next likely illustrate the relative performance, even without performing a
much more careful tuning of the code. The timing results shown here were obtained in a
2.9 Ghz Intel Core i7 processor for the two-dimensional case, and 6-core Intel Westmere
in the three-dimensional case. The same base code with appropriate modifications was
used to implement the three methods, and no special efforts were made to optimize any
particular part of the code. In all cases, no preconditioner and GMRES were used to
solve each system of equations.
The results are shown in Figs. 23 and 24. The timing results for HDG are clearly
encouraging, especially for the computation of the approximation of the gradient of the
solution, and hence the stresses. The two curves cross just before the second subdivision
in three-dimensions which, again, are only a few tens of thousands degrees of freedom.
6. Summary and Conclusions
We have shown how to recast the method proposed in [2, 1] as a minimization problem
and that if the stabilization function is not big enough, convergence to the wrong solution
might take place. We have also shown numerical experiments displaying the performance
of the method in several trying situations in two- and, for the first time, three-space
dimensions and concluded that the method behaves fairly well, especially compared with
the CG and the DG methods.
28
1
1
10
10
0
0
10
10
−1
−1
Error
Error
10
−2
10
10
−2
10
−3
−3
10
10
HDG
HDG
CG
DG
CG
DG
−4
10
1
−4
2
10
10
3
10
4
10
5
10
10
6
10
1
2
10
10
3
10
4
10
5
10
6
10
7
10
Non-zero entries in matrix
Size of the system
b h )||L2 (B0 ) versus computational cost for a two-dimensional geometry
Figure 20: ||∇ϕ − DDG (ϕh , ϕ
−1
−1
10
10
HDG
HDG
CG
CG
−2
−2
10
Error
Error
10
−3
−3
10
10
−4
10
−4
2
3
10
10
4
10
5
10
6
10
10
7
10
3
4
10
10
5
10
6
10
7
10
8
10
9
10
Non-zero entries in matrix
Size of the system
Figure 21: ||ϕ − ϕh ||L2 (B0 ) versus computational cost for a three-dimensional geometry.
0
0
10
10
HDG
HDG
CG
CG
−1
−1
10
Error
Error
10
−2
−2
10
10
−3
10
−3
2
10
3
10
4
10
5
10
6
10
10
7
10
3
10
4
10
5
10
6
10
7
10
8
10
9
10
Non-zero entries in matrix
Size of the system
b h )||L2 (B0 ) versus computational cost for a three-dimensional geometry.
Figure 22: ||∇ϕ − DDG (ϕh , ϕ
29
1
0
10
10
HDG
CG
DG
0
−1
10
−2
10
10
−1
10
Error
Error
HDG
CG
DG
−3
−2
10
10
−3
−4
10
10
−4
−5
10 −3
10
−2
10
−1
10
0
10
1
10
2
10
10 −3
10
3
10
−2
10
−1
10
0
10
1
10
2
10
3
10
Time
Time
b h )||L2 (B0 )
(b) ||∇ϕ − DDG (ϕh , ϕ
(a) ||ϕ − ϕh ||L2 (B0 )
Figure 23: Error of the solution versus time taken to solve the non-linear system of equations for a
two-dimensional geometry.
0
−1
10
10
HDG
HDG
CG
CG
−1
−2
10
Error
Error
10
−2
−3
10
10
−3
−4
10 −2
10
−1
10
0
10
1
10
2
10
3
10
4
10
10 −2
10
−1
10
0
10
1
10
2
10
3
10
4
10
Time
Time
b h )||L2 (B0 )
(b) ||∇ϕ − DDG (ϕh , ϕ
(a) ||ϕ − ϕh ||L2 (B0 )
Figure 24: Error of the solution versus time taken to solve the non-linear system of equations for a
three-dimensional geometry.
30
Let us end by pointing out that the characterization of the HDG method [2, 1] as
a minimization problem can be easily extended to many other problems of continuum
mechanics that can also be formulated as minimization problems. How to take advantage
of this fact, and in particular, how to use it to define the stabilization function τ in an
automatic manner, constitutes the subject of ongoing work.
Acknowledgements
This work was supported by the Franklin P. Johnson Jr. Stanford Graduate Fellowship to Hardik Kabaria. Adrian J. Lew is grateful for the support for this work from the
following sources: Department of the Army Research Grant, grant number: W911NF07- 2-0027; NSF Career Award, grant number: CMMI-0747089; and NSF, grant number
CMMI-1301396. Bernardo Cockburn would like to acknowledge the support of the NSF
through the DMS-1115331 grant.
References
[1] S.-C. Soon, Hybridizable discontinuous Galerkin methods for solid mechanics, Ph.D. thesis, University of Minnesota (2008).
[2] S.-C. Soon, B. Cockburn, H. K. Stolarski, A hybridizable discontinuous Galerkin method for linear
elasticity, Internat. J. Numer. Methods Engrg. 80 (8) (2009) 1058–1092. doi:10.1002/nme.2646.
URL http://dx.doi.org/10.1002/nme.2646
[3] P. Lew, A.and Neff, D. Sulsky, M. Ortiz, Optimal BV estimates for a discontinuous Galerkin method
for linear elasticity, AMRX Appl. Math. Res. Express (3) (2004) 73–106.
[4] P. Hansbo, M. Larson, Discontinuous finite element methods for incompressible and nearly incompressible elasticity by use of Nitsche’s method, Comput. Methods Appl. Mech. Engrg. 191 (2002)
1895–1908.
[5] B. Cockburn, D. Sch¨
otzau, J. Wang, Discontinuous Galerkin methods for incompressible elastic
materials, Comput. Methods Appl. Mech. Engrg. 195 (2006) 3184–3204, c. Dawson, Ed.
[6] F. Celiker, B. Cockburn, S. G¨
uzey, R. Kannapady, S.-C. Soon, H. K. Stolarski, K. K. Tamma,
Discontinuous Galerkin methods for Timoshenko beams, in: Numerical Mathematics and Advanced
Applications, ENUMATH 2003, Springer, 2003, pp. 221–231.
[7] F. Celiker, B. Cockburn, H. Stolarski, Locking-free optimal discontinuous Galerkin methods for
Timoshenko beams, SIAM J. Numer. Anal. 44 (2006) 2297–2325.
[8] G. Engel, K. Garikipati, J. T. R. Hughes, M. Larson, L. Mazzei, R. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and
continuum mechanics with applications to thin beams and plates, and strain gradient elasticity,
Comput. Methods Appl. Mech. Engrg. 191 (2002) 3669–3750.
[9] S. Brenner, L.-Y. Sung, C 0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, J. Sci. Comput. 22/23 (2005) 83–118.
[10] D. N. Arnold, F. Brezzi, L. D. Marini, A family of discontinuous Galerkin finite elements for the
Reissner-Mindlin plate, J. Sci. Comput. 22/23 (2005) 25–45.
[11] S. G¨
uzey, B. Cockburn, H. Stolarski, The embedded discontinuous Galerkin methods: Application
to linear shells problems, Internat. J. Numer. Methods Engrg. 70 (2007) 757–790.
[12] A. Ten Eyck, A. Lew, Discontinuous Galerkin methods for non-linear elasticity, Internat. J. Numer.
Methods Engrg. 67 (2006) 1204–1243.
[13] L. Noels, R. Radovitzky, A general discontinuous Galerkin method for finite hyperelasticity. Formulation and numerical applications, Internat. J. Numer. Methods Engrg. 68 (1) (2006) 64–97.
[14] B. Cockburn, C.-W. Shu (Eds.), Special issue on discontinuous Galerkin methods, Vol. 22 and 23
of J. Sci. Comput., Springer, 2005.
[15] C. Dawson (Ed.), Special issue on discontinuous Galerkin methods, Vol. 195 of Comput. Methods
Appl. Mech. Engrg., Elsevier, 2006.
[16] R. Gracie, H. Wang, T. Belytschko, Blending in the extended finite element method by discontinuous
Galerkin and assumed strain methods, Internat. J. Numer. Methods Engrg. 74 (11) (2008) 1645–
1669.
31
[17] Y. Shen, A. Lew, An optimally convergent discontinuous galerkin-based extended finite element
method for fracture mechanics, International journal for numerical methods in engineering 82 (6)
(2010) 716–755.
[18] Y. Shen, A. Lew, Stability and convergence proofs for a discontinuous-galerkin-based extended finite
element method for fracture mechanics, Computer Methods in Applied Mechanics and Engineering
199 (37) (2010) 2360–2382.
[19] J. K. Djoko, F. Ebobisse, A. T. McBride, B. D. Reddy, A discontinuous Galerkin formulation for
classical and gradient plasticity. I. Formulation and analysis, Comput. Methods Appl. Mech. Engrg.
196 (2007) 3881–3897.
[20] J. K. Djoko, F. Ebobisse, A. T. McBride, B. D. Reddy, A discontinuous Galerkin formulation for
classical and gradient plasticity. II. Algorithms and numerical analysis, Comput. Methods Appl.
Mech. Engrg. 197 (1-4) (2007) 1–21.
[21] A. T. Eyck, F. Celiker, A. Lew, Adaptive stabilization of discontinuous Galerkin methods for
nonlinear elasticity: Motivation, formulation, and numerical examples, Comput. Methods Appl.
Mech. Engrg. 197 (2007) 1–21.
[22] A. Eyck, F. Celiker, A. Lew, Adaptive stabilization of discontinuous Galerkin methods for nonlinear
elasticity: analytical estimates, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2989–3000.
[23] L. Noels, R. Radovitzky, Alternative approaches for the derivation of discontinuous Galerkin methods for nonlinear mechanics, Journal of Applied Mechanics 74 (2007) 1031–1036.
[24] L. Noels, R. Radovitzky, A new discontinuous Galerkin method for kirchhoff-love shells, Internat.
J. Numer. Methods Engrg. 197 (2008) 2901–2929.
[25] L. Noels, R. Radovitzky, An explicit discontinuous Galerkin method for non-linear solid dynamics: Formulation, parallel implementation and scalability properties, Internat. J. Numer. Methods
Engrg. 74 (2008) 1393–1420.
[26] G. Fu, B. Cockburn, K. Shi, Analysis of an HDG method for linear elasticitySubmitted.
[27] N. Nguyen, J. Peraire, B. Cockburn, High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics, J. Comput. Phys. 230 (2011) 3695–3718.
[28] N. C. Nguyen, J. Peraire, Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics., J. Comput. Physics 231 (18) (2012) 5955–5988.
[29] B. Cockburn, J. Gopalakrishnan, N. Nguyen, J. Peraire, F. Sayas, Analysis of an HDG method for
Stokes flow, Math. Comp. 80 (2011) 723–760.
[30] N. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous Galerkin
method for the incompressible Navier-Stokes equations, J. Comput. Phys. 230 (2011) 1147–1170.
[31] B. Cockburn, K. Shi, Superconvergent HDG methods for linear elasticity with weakly symmetric
stresses, IMA J. Numer. Anal. 33 (2013) 747–770.
[32] A. Buffa, C. Ortner, Compact embeddings of broken Sobolev spaces and applications, IMA J. Num.
Anal. 29 (2009) 827–855.
[33] J. Marsden, T. Hughes, Mathematical foundations of elasticity, Dover Publications, 1994.
[34] J. Ball, Some open problems in elasticity, Geometry, mechanics, and dynamics (2002) 3–59.
[35] D. N. Arnold, F. Brezzi, B. Cockburn, L. D. Marini, Unified analysis of discontinuous Galerkin
methods for elliptic problems, SIAM journal on numerical analysis 39 (5) (2002) 1749–1779.
[36] B. Dacorogna, Direct methods in the calculus of variations, Vol. 78, Springer, 2008.
[37] C. Mardare, A nonlinear Korn inequality with boundary conditions and its relation to the existence
of minimizers in nonlinear elasticity, Comptes Rendus Mathematique 349 (34) (2011) 229 – 232.
doi:http://dx.doi.org/10.1016/j.crma.2011.01.011.
[38] A. Gent, P. Lindley, Internal rupture of bonded rubber cylinders in tension, Proceedings of the
Royal Society of London. Series A. Mathematical and Physical Sciences 249 (1257) (1959) 195–205.
[39] X. Xu, D. Henao, An efficient numerical method for cavitation in nonlinear elasticity, Mathematical
Models and Methods in Applied Sciences 21 (08) (2011) 1733–1760.
[40] Y. Lian, Z. Li, A numerical study on cavitation in nonlinear elasticitydefects and configurational
forces, Mathematical Models and Methods in Applied Sciences 21 (12) (2011) 2551–2574.
[41] T. Nakamura, O. Lopez-Pamies, A finite element approach to study cavitation instabilities in nonlinear elastic solids under general loading conditions, International Journal of Non-Linear Mechanics
47 (2) (2012) 331–340.
[42] D. Ruprecht, H. M¨
uller, A scheme for edge-based adaptive tetrahedron subdivision, in: H.-C. Hege,
K. Polthier (Eds.), Visualization and Mathematics, Springer Verlag, Heidelberg, 1998, pp. 61–70.
32
A. Solving the Non-Linear Equations
We describe next our implementation of the method. To solve the system of non-linear
equations that result from the HDG discretization, we adopted a standard Newton’s
method. Alternative algorithms that take advantage of the structure of the problem as a
minimization were not explored here, but could be promising. In particular, a non-linear
conjugate gradients algorithm does not need the computation of the second variations
of Ih , namely, the stiffness matrix. Finally, the potential energy (2) may have multiple
minimizers, so often a version of a continuation method is needed to precisely define the
minimizers sought. Such methods were not used here.
A.1. Newton’s Method
b h ) of the functional Ih [ϕh , ϕ
b h]
We know, by Theorem 1, that the critical points (ϕh , ϕ
b h (g) must satisfy the following equation:
over Vh × V
∂W
h
h
b )), DDG (v, µ)
b h , v − µi∂Th = (f , v)Th + hT, µi∂N B0
(DDG (ϕ , ϕ
+ τ h ϕh − ϕ
∂F
Th
b h (0)
for all (v, µ) ∈ Vh × V
To apply Newton’s method to this system of equations, it is convenient to obtain an
explicit expression for the DG-derivative.
b h ). It is convenient to introduce the lifting operator
Explicit computation of DDG (ϕh , ϕ
R : [L2 (∂Th )]d → Gh (see, e.g., [35, 12]), which satisfies that
(R(b
v), G)K = hb
v, Gni∂K
∀ G ∈ Gh .
(24)
It then follows that, after using integration by parts, (6) can be restated as
b h ), G)Th = (∇ϕh , G)Th + hb
(DDG (ϕh , ϕ
ϕh − ϕh , Gni∂Th
= (∇ϕh + R(b
ϕh − ϕh ), G)Th
for all ∀ G ∈ Gh , from where we obtain an explicit formula for the DG-derivative,
namely,
b h ) = Dh = ∇ϕh + R(b
DDG (ϕh , ϕ
ϕh − ϕh ).
(25)
The advantage of introducing the lifting operator is that it is possible to compute
b h , since it is a
basis functions for it that are conjugate to the degrees of freedom in V
linear operator. In this way, it is not necessary to solve the system (6) every time a DGderivative is needed. These basis functions can then be used to construct R(b
ϕh − ϕh )
h
ba ∈ V
b h , the
b − ϕh . More specifically, given a basis function N
for any value of ϕ
corresponding basis function for the lifting operator Ra ∈ Gh is defined by
b a , Gni∂K
(Ra , G)K = hN
∀G ∈ Gh .
(26)
This is an element-by-element computation, which involves inverting the elemental mass
matrix for Gh . These basis functions can be computed for any given mesh, and used for
all iterations
in the non-linear solver, and for all problems, in such mesh. Then, given
P
ba ∈ V
b h , the lifting operator follows as
µ = a µa N
X
R(µ) =
µa Ra .
(27)
a
33
b h` }`≥0 be the sequence generated by Newton’s method.
Newton’s method step. Let {ϕh` , ϕ
h
b ` ) := (ϕh`+1 , ϕ
b h`+1 ) − (ϕh` , ϕ
b h` ) is defined by the linearized
Then, the increment (δϕh` , δ ϕ
equation
∂2W
b h` )) : DDG (δϕh` , δ ϕ
b h` )), DDG (v, µ)
(DDG (ϕh` , ϕ
∂2F
+
Th
b h` , v − µi∂Th = −R` (v, µ)
τ h δϕh` − δ ϕ
(28)
b h (0), where
for all (v, µ) ∈ Vh × V
R` (v, µ) :=
∂W
b h` )), DDG (v, µ)
(DDG (ϕh` , ϕ
∂F
Th
b h` , v − µi∂Th − (f , v)Th − hT, µi∂N B0 .
+ τ h ϕh` − ϕ
Equation (28) defines a linear system of equations to be solved at each Newton’s step.
The matrix of such system is the stiffness matrix of the problem, and it changes at each
Newton’s step. To solve this system, it is again possible to exploit the ability of the
HDG method to statically condense the degrees of freedom for δϕh` inside each element,
b h` only. The value of δϕh` is recovered later.
to solve the global system of equations for δ ϕ
This feature can be clearly seen from the sparsity pattern of the stiffness matrix, see Fig.
1b. As shown there, the submatrix corresponding to the degrees of freedom for δϕh` is
block diagonal, with a dense diagonal block per element. The coupling of these degrees of
b h` , and hence, as mentioned, can be statically
freedom only occurs through the fluxes δ ϕ
condensed. We did not implement this approach here, but still obtained very encouraging
results, as discussed later. For a more detailed discussion of this implementation, see
[1, 2, 28].
B. Proof of Lemma 1
proof. A straightforward computation gives us that
∂
∂
b h, ]
b h )),
b h, )
= ∂W/∂F(DDG (ϕh , ϕ
Ih [ϕh, , ϕ
DDG (ϕh, , ϕ
∂
∂
=0
=0 Th
(29)
b h i∂N B0
− (f , δϕh )Th − hT, δ ϕ
b h ), (ϕh − ϕ
b h )i∂Th .
+ h τ (δϕh − δ ϕ
To calculate the first term of the right-hand side, we only have to consider the definition
of the DG-derivative (6) to obtain
∂
h, h,
b )
b h , Gni∂Th
DDG (ϕ , ϕ
,G
= −(δϕh , ∇ · G)Th + hδ ϕ
∀ G ∈ Gh .
∂
=0
Th
b h ), G
= DDG (δϕh , δ ϕ
Th
34
It then follows that
∂
h, h,
b )
b h)
= DDG (δϕh , δ ϕ
DDG (ϕ , ϕ
∂
=0
(30)
b h (0), and the first equality of the lemma is obtained after
b h ) ∈ Vh × V
for all (δϕh , δ ϕ
replacing in (29).
For the second equality, we use the definition of Ph and (30) to obtain
∂
b h )),
b h, )
DDG (ϕh, , ϕ
Θ := ∂W/∂F(DDG (ϕh , ϕ
∂
=0 Th
h
h
h
b )
= P , DDG (δϕ , δ ϕ
Th
b h , Ph ni∂Th ,
= −(δϕh , ∇ · Ph )Th + hδ ϕ
b h − δϕh i∂Th ,
= (Ph , ∇δϕh )Th + hPh n, δ ϕ
after a simple integration by parts. Then, we can conclude that
∂
h, h,
b ]
b h − δϕh i∂Th
= (Ph , ∇δϕh )Th + hPh n, δ ϕ
Ih [ϕ , ϕ
∂
=0
b h i∂N B0
− (f , δϕh )Th − hT, δ ϕ
b h ), (ϕh − ϕ
b h )i∂Th ,
+ h τ (δϕh − δ ϕ
b h n and
and, after using the definition of the normal component of the numerical trace P
rearranging terms, that
∂
b h n, δϕh − δ ϕ
b h, ]
b h i∂Th ,
b h i∂N B0 −h P
= (Ph , ∇δϕh )Th −(f , δϕh )Th −hT, δ ϕ
Ih [ϕh, , ϕ
∂
=0
and the second equality follows. This completes the proof of Lemma 1.
35

An HDG Method for unsteady compressible flows - IGPM

Axone Single-Axis Tracker Cut Sheet

Multi-physics imaging models in discontinuous media

Aluminium DVA mesh,Supa mesh,limited vision Mesh

Solaris-Shield

numerical modeling of complex 3d woven composite

Managerial Economics - Dronacharya College of Engineering

cal-tb133 chart

JRF Position in a DST Project

Answers to Webinar "Wind farm flow modelling using CFD