Because of their relevance in the development of new technologies bending theories for thin sheets have attracted considerable attention within applied mathematics in the last two decades, starting from the seminal article [FrJaMu02]. In this article the folding of thin elastic sheets along a prepared curved arc is considered which naturally leads to bending effects, cf. Figure 1. This setting has only been partially addressed mathematically, e.g., [DunDun82], but has recently attracted considerable attention in applied sciences, cf., e.g. [Specketal15, CheMah19, ChDuMa19, PeHaLa19-book, CGGP19, LiPlJaFe21] and references therein. Particular applications arise in the design of cardboxes and bistable switching devices that make use of corresponding flapping mechanisms. It is our aim to derive a mathematical description via a rigorous dimension reduction from three-dimensional hyperelasticity and to devise effective numerical methods that correctly predict large deformations in practically relevant settings.
To describe our approach we let be a bounded Lipschitz domain that represents the sheet without thickness and be a curve that models the crease, i.e., the arc along which the sheet is folded. The corresponding three-dimensional model involves a thickness parameter and the thin domain . The material is weakened (or damaged) in a neighbourhood of width around the arc and we consider for given functions and the hyperelastic energy functional for a deformation
Here is small with value close to the arc and approximately away from the arc, the function is a typical free energy density. Hence, the factor models a reduced elastic response of the material close to . By appropriately relating the thickness , the intactness fraction , and width of the prepared region, we obtain for a meaningful dimensionally reduced model which seeks a minimizing deformation for the functional
where is the second fundamental form related to the parametrization which is weakly differentiable in with second weak derivatives away from , i.e., we consider
Moreover, is required to satisfy the pointwise isometry condition
which implies that no shearing and stretching effects occur. The minimization of is supplemented by boundary conditions and possible body forces. Note that for our scaling of the parameters , , and , e.g., and , no energy contributions such as a penalization of the folding angle arise from the crease, i.e., the material can freely fold along the arc and is in general only continuous on .
In the absence of a crease the model coincides with Kirchhoff’s plate bending functional, e.g., describing the deformation of paper, rigorously derived in [FrJaMu02]. We slightly modify the arguments given there to take account of the fact that any asymptotic deformation which does not only bend but also folds has infinite Kirchhoff energy. This is because the Hessian corresponding to a folding deformation is not square integrable. Other variants of the setting from [FrJaMu02] have been addressed, in [FJMM03, FrJaMu02b, HorVel18, Velc15, Schm07] and many others.
A typical setting and an experiment are shown in Figure 1. Interesting phenomena phenomena take place at the folding arc. It was observed in [DunDun82] that a deformation on one side of the arc locally restricts possible deformations on the other side. In fact, only a finite number of scenarios is possible and either the gradients coincide along the arc resulting in a smooth deformation or a discontinuity occurs and the deformation is locally up to a sign uniquely determined. This important effect is a result of the isometry condition and the related physical property that thin elastic sheets are unshearable in the bending regime. Moreover, this effect arises in biology and inspires the development of new technologies and design in architecture [Specketal15].
Our numerical method to approximate minimizers is based on a discontinuous Galerkin method from [BoNoNt21] which generalizes the approach based on a nonconforming method of [Bart13a]. It allows us to define a discrete curve as a union of element sides that approximates and to account for possible discontinuities of deformation gradients along by simply removing typical jump terms in the discrete formulation. The continiuity condition on the deformation along the interface is similar to a simple support boundary condition in linear bending theories. For such problems the plate paradox, cf. [BabPit90], states that convergence of approximations resulting from polyhedral approximations of curved domains or interfaces fails in general. We therefore consider piecewise quadratic approximations of . We note however that we do not observe significant differences to approximations obtained with piecewise linear arcs in the nonlinear setting under consideration.
We assume for simplicity and without loss of generality that and use that the Frobenius norm of the second fundamental form equals the Frobenius norm of the Hessian in the case of an isometry , i.e., . Our numerical method then uses discontinuous deformations from an isoparametric finite element space subordinated to a partitioning of and a suitably defined discrete Hessian in which define the discrete energy functional
The last two terms are typical stabilization terms with the mesh-size function on the skeleton of the triangulation that guarantee coercivity and enforce continuity of and the elementwise gradient across interelement sides or essential boundary conditions on certain boundary sides as the mesh-size tends to zero. The union of all such sides is the set of active sides which is denoted by . Crucial here is that the penalized continuity of is not imposed along the discrete folding arc and that related consistency terms do not enter in the definition of the discrete Hessian . The important isometry condition is imposed up to a tolerance via averages on elements, i.e., we require that
for a suitable tolerance . To obtain accurate approximations of large deformations choosing a small parameter is desirable. Its particular choice is dictated by available density results. If, e.g., the density of smooth folding isometries can be guaranteed similar to [Horn11b] even a pointwise controlled violation criterion can be used; otherwise a condition has to be satisfied, cf. [BoNoNt21]. Boundary conditions on a part of the boundary are included in the discrete problem via an appropriate definition of the jump terms. A justification of the discrete energy functionals via convergence as in [Bart13a, BaBoNo17, BoNoNt21] is in preparation.
The iterative solution of the constrained minimization problem follows the ideas of [Bart13a, BaBoNo17, BoNoNt21] and is realized by a discrete gradient flow with a suitable linearization of the constraint. In particular, we consider the linear space of variations for a given deformation defined via
We furthermore let be an inner product and a step size. Given an approximation we look for a correction such that
for all and with the discrete bilinear form associated with the discrete quadratic energy functional . The existence of a unique solution is an immediate consequence of the Lax–Milgram lemma and we define the new approximation
This implies the interpretation of the symbol as a backward difference quotient. By choosing one directly obtains the energy decay property for ,
in particular, we see that as , i.e., the iteration becomes stationary. Because of the orthogonality relation included in the space we can bound the violation of the isometry constraint by repeatedly replacing , i.e., for all we have
If the discrete gradient flow metric controls the
norm of the elementwise gradient, we obtain from the energy decay law the estimate
where is the initial isometry violation. In particular, we see that if and are sufficiently small, an arbitrary accuracy can be achieved and we have independently of the number of iterations .
Using the numerical scheme we simulate various scenarios that are motivated by practical applications, e.g., how the shape of the folding arc affects the flapping mechanism, or address subtle analytical features of solutions, e.g., the occurrence of energy concentrations when the curve has a kink. Besides that we illustrate the robustness of the numerical method with respect to the choice of stabilization parameters and discuss the construction of suitable deformations that serve as starting values in the discrete gradient flow. Our experiments show that large deformations in highly nontrivial settings can be accurately computed with moderate resolution.
The outline of the article is as follows. In Section 2 we describe the general setup and the dimensionally reduced model. Its rigorous derivation is given in Section 3. The discontinuous Galerkin finite element method is derived and stated in Section 4. Numerical experiments are reported in Section 5.
2.1. Hyperelasticity for plates
For a bounded Lipschitz domain we consider a plate of thickness occupying the domain in the reference configuration. The elastic energy stored in the configuration determined by a deformation is given by
Here is a frame indifferent stored energy function and as in [FrJaMu02] we impose the following conditions:
and in a neighbourhood of .
is frame indifferent, i.e., for all and all . Moreover, .
There is a constant such that .
To analyze the limiting behaviour as it is convenient to work on the fixed domain
where . We define a rescaled deformation by setting . Then the (re-scaled) elastic energy is given by
where and .
Throughout this article we use standard notation related to Sobolev spaces, e.g., denotes the set of times weakly differentiable, -valued functions in whose weak derivatives are -integrable. norms are often used without specifying a domain when there is no ambiguity, and we abbreviate the norm on by . We occasionally omit target domains when this is clear from the context. The open ball of radius around a point is denoted by . For integral functionals ocurring below, it is often useful to specify their integration domains explicitly, e.g., we write
For the canonical choice, e.g., , this argument is usually omitted. For we identify maps defined on with their trivial extension to .
2.3. Folded plates
Our aim is to modify the arguments from [FrJaMu02] in order to allow for folding effects along a prescribed curve, see Figure 1. In applications, the folding curve is prescribed by weakening the material along it. For a mathematical description we let be a smooth embedded Lipschitz curve intersecting the boundary at both endpoints. Such a curve disconnects into precisely two connected components and of . We assume, moreover, that intersects the boundary of transversally and that it is such that both and become Lipschitz domains. In particular, is locally a Lipschitz graph. For we define
As explained in the introduction, we let , be parameters that define the width of the prepared region and the amount of the material intactness. We then define by
denotes the characteristic function of a set. With this we consider the (re-scaled) three-dimensional energy functional
Our limit passage for leads to a pointwise isometry constraint and discontinuities of gradients across the arc , i.e., we are led to the set of asymptotically admissible deformations
The corresponding asymptotic energy functional is defined as
Here, is the second fundamental form of the surface parametrized by with unit normal , i.e.,
and is obtained by relaxing, over the third column and row, the quadratic form corresponding to the Hessian of
at the identity matrix, i.e.,
where for given the matrix is obtained by consistently appending a row and column defined by . Note that by hypotheses (H1)-(H3) we have the Taylor expansion
for and .
(i) Observe that every belongs to
, because and its gradient
is bounded almost everywhere on . In particular, is continuous
(ii) Notice that a Lipschitz function is in precisely if there is an such that
The purpose of this section is to prove the following result.
Let be a bounded Lipschitz domain, let be a smooth embedded curve with both endpoints on intersecting transversally at its endpoints. Let , be such that
as well as
Assume that are such that
Then there exists a subsequence (not relabelled) and such that weakly in and locally strongly in .
Assume that weakly in . Then
Let . Then there exist such that
This theorem is proven in the next two sections. It relies on arguments and results in [FrJaMu02].
3.1. Compactness and lower bound
In what follows we will frequently omit the index in and .
We omit the index in and . By the definition of in (3) we have
On the other hand, by (6),
We conclude that .
Hence is uniformly bounded in due to (2.1),
so there exists such that weakly
in , after taking subsequences.
Since and , we have
Now fix and define for , ; both are Lipschitz domains. We can apply [FrJaMu02, Theorems 4.1 and 6.1 (i)] on each . Hence strongly in and with
Here is the second fundamental form of on . Since is an isometric immersion, by [FJM2, Proposition 6] we have
Hence (9) implies that
where does not depend on . Since is locally a Lipschitz graph, the area of converges to as . Hence and
Summarising, we have
locally strongly in
According to (9), for every we have
The right-hand side does not depend on . Taking the supremum over all small , we therefore see that
which proves the proposition. ∎
3.2. Recovery sequence
In the proof of the upper bound in Theorem 3.1 we will use the following lemma.
Let be a bounded Lipschitz domain. Then there exists a constant , depending only on the Lipschitz constant of , such that the following is true: if satisfies and if we set
then intersects for each .
This is standard; we include the proof
By Definition 1.3 in [Giaquinta, Chapter III] and the remark
following it, there exists a constant , depending only
on the Lipschitz constant of , such that
whenever and .
If did not intersect , then and thus we would have
a contradiction. ∎
The main result of this section is the following proposition.
As before denote the connected components of .
Denote the restriction of to by .
In this proof we will omit the index and write instead of
and instead of . The letter denotes constants that do
not depend on as .
For each let be a smooth cutoff function with on and on ; we choose it such that
Set and . Denote by the normal to and define
Recall that each is a Lipschitz domain. So
we can extend each to a map in
We extend to a map in .
As in the proof of [FrJaMu02, Theorem 6.1 (ii)] we truncate these maps and thus obtain maps and satisfying the bound
while at the same time there is a set with
Let and define the recovery sequence
By (13) we have
Recalling that and , we see that on
We claim that
In fact, let . By definition of there exists such that . Since on and since is Lipschitz, we have
Since and since is Lipschitz, we have
After repeating the argument with and , we conclude that
Therefore, since due to the regularity of ,
The right-hand side converges to zero due to (8).