 # Modeling and simulation of thin sheet folding

The article addresses the mathematical modeling of the folding of a thin elastic sheet along a prescribed curved arc. A rigorous model reduction from a general hyperelastic material description is carried out under appropriate scaling conditions on the energy and the geometric properties of the folding arc in dependence on the small sheet thickness. The resulting two-dimensional model is a piecewise nonlinear Kirchhoff plate bending model with a continuity condition at the folding arc. A discontinuous Galerkin method and an iterative scheme are devised for the accurate numerical approximation of large deformations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

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

 Eh(z)=∫Ωhfε,r(x)W(∇z)dx.

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

 EK(y)=124∫S∖ΣQ(A)dx′,

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

 y∈W2,2(S∖Σ;R3)∩W1,∞(S;R3).

Moreover, is required to satisfy the pointwise isometry condition

 (∇y)T(∇y)=I,

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. Figure 1. Left: Geometry of a prepared elastic sheet S with folding arc Σ that separates regions S1 andf S2; boundary conditions are imposed in the points xB,x′B. Right: When the boundary points xB and x′B are moved towards each other a flapping mechanism occurs.

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

 ˜EK(˜Y)=124∫S∖Σ|˜H(˜Y)|2dx′+γ02∫∪˜Eah−3˜E|⟦˜Y⟧|2ds+γ12∫∪˜Ea∖˜Σh−1˜E|⟦˜∇˜Y⟧|2ds.

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

 (1) ˜Y∈˜A={˜Z∈˜V3:∑T∈˜T∣∣∫T(˜∇˜Z)T(˜∇˜Z)−Idx′∣∣≤˜ϱ}

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

 ˜F[˜Z]={˜W∈˜V3:∫T(˜∇˜Z)T(˜∇˜W)+(˜∇˜W)T(˜∇˜Z)dx′=0 for all T∈˜T}.

We furthermore let be an inner product and a step size. Given an approximation we look for a correction such that

 (2) (dt˜Yk,˜V)⋆+˜aK(˜Yk−1+τdt˜Yk,˜V)=0

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

 ˜Yk=˜Yk−1+τdt˜Yk−1.

This implies the interpretation of the symbol as a backward difference quotient. By choosing one directly obtains the energy decay property for ,

 ˜EK[˜Yℓ]+τℓ∑k=1∥dt˜Yk∥2⋆≤˜EK[˜Y0]=˜E0K

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

 ∫T[˜∇˜Yℓ]T˜∇˜Yℓ−Idx′=∫T[˜∇˜Y0]T˜∇˜Y0−I+τ2ℓ∑k=1[˜∇dt˜Yk]T˜∇dt˜Ykdx′.

If the discrete gradient flow metric controls the

norm of the elementwise gradient, we obtain from the energy decay law the estimate

 ∑T∈˜T∣∣∫T[˜∇˜Yℓ]T˜∇˜Yℓ−Idx′∣∣≤ε0+c⋆˜E0Kτ,

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. Preliminaries

### 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

 ∫ΩhW(∇z)dx.

Here is a frame indifferent stored energy function and as in [FrJaMu02] we impose the following conditions:

1. and in a neighbourhood of .

2. is frame indifferent, i.e., for all and all . Moreover, .

3. There is a constant such that .

To analyze the limiting behaviour as it is convenient to work on the fixed domain

 Ω=S×I,

where . We define a rescaled deformation by setting . Then the (re-scaled) elastic energy is given by

 ˜Eh(yh)=∫ΩW(∇hyh)dx,

where and .

### 2.2. Notation

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

 E(y;~S)=∫~SF(y)dx′.

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

 ΣR=BR(Σ)=⋃x∈ΣBR(x).

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

 (3) fh=εhχΣrh+1−χΣrh,

where

denotes the characteristic function of a set

. With this we consider the (re-scaled) three-dimensional energy functional

 (4) Eh(y)=∫Ωfh(x′)W(∇hyh(x))dx.

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

 A(S,Σ)={u∈W1,2(S;R3)∩W2,2(S∖Σ;R3):(∇u)T(∇u)=I a.e. on S}.

The corresponding asymptotic energy functional is defined as

 (5) EK(y)={124∫S∖ΣQ(A)dx′ if u∈A(S,Σ),+∞ otherwise.

Here, is the second fundamental form of the surface parametrized by with unit normal , i.e.,

 A=(∇n)T(∇y),

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.,

 Q(A)=mind∈R3D2W(I)[(A|d),(A|d)],

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

 W(I+hF)=h22D2W(I)[F,F]+o(h2)

for and .

###### Remarks 2.1.

(i) Observe that every belongs to , because and its gradient is bounded almost everywhere on . In particular, is continuous on .
(ii) Notice that a Lipschitz function is in precisely if there is an such that

 ∇2f=F in the sense of distributions on S∖Σ.

## 3. Gamma-Convergence

The purpose of this section is to prove the following result.

###### Theorem 3.1.

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

 (6) limsuph→0h2εh<∞

and

 (7) limsuph→0hrh<∞

as well as

 (8) limsuph→0εhrhh2=0.

Define as in (3), define as in (4) and define as in (5).
Then deformations with finite bending energy are compact and Gamma-converges to . More precisely:

1. Assume that are such that

 limsuph→01h2Eh(yh)<∞.

Then there exists a subsequence (not relabelled) and such that weakly in and locally strongly in .

2. Assume that weakly in . Then

 EK(y)≤liminfh→0Eh(yh).
3. Let . Then there exist such that

 limh→01h2Eh(yh)=EK(y).

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 .

###### Proposition 3.2.

Let and as in Theorem 3.1, let as and let satisfy (6). If satisfy

 limsuph→01h2Eh(yh)<∞,

then there exists an isometric immersion

 y∈W1,2(S,R3)∩W2,2(S∖Σ,R3)

such that, after taking subsequences, weakly in and locally strongly in as . Moreover,

 EK(y;Si)≤liminfh→01h2Eh(yh;Si×I) for i=1,2.
###### Proof.

We omit the index in and . By the definition of in (3) we have

 ∫Ω∖(Σr×I)W(∇hyh)dx≤Ch2≤C.

On the other hand, by (6),

 ε∫Σr×IW(∇hyh)dx=∫Σr×IfhW(∇hyh)dx≤Ch2≤Cε.

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

 limsupR→0∫Ω∖(ΣR×I)W(∇hyh)dx≤Ch2.

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

 (9) 124∫SRiQ(A)dx′≤liminfh→0h−2∫SRi×IW(∇hyh)dx≤C.

Here is the second fundamental form of on . Since is an isometric immersion, by [FJM2, Proposition 6] we have

 |∇2y|=|A| almost everywhere on S∖ΣR.

Hence (9) implies that

 (10) ∥∇2y∥L2(SR1)+∥∇2y∥L2(SR2)≤C,

where does not depend on . Since is locally a Lipschitz graph, the area of converges to as . Hence and

 ∫SiQ(A)dx′=limsupR→0∫SRiQ(A)dx′<∞.

Summarising, we have locally strongly in and .
According to (9), for every we have

 EK(y;SRi) ≤liminfh→0h−2Eh(yh;SRi×I) ≤liminfh→0h−2Eh(yh;Si×I).

The right-hand side does not depend on . Taking the supremum over all small , we therefore see that

 EK(y;Si)≤liminfh→0h−2Eh(yh;Si×I),

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.

###### Lemma 3.3.

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

 R=√2|M|δ,

then intersects for each .

###### Proof.

This is standard; we include the proof for convenience. 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

 2|M|=δR2≤|BR(x)∩U|≤|M|,

The main result of this section is the following proposition.

###### Proposition 3.4.

Let and as in Theorem 3.1, let , satisfy (7) and (8) and let . Then there exist such that weakly in and as .

###### Proof.

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

 (11) ∥∇′ηh∥L∞(Σr)≤Cr.

Set and . Denote by the normal to and define

 nh=ηh1n1+ηh2n2.

Recall that each is a Lipschitz domain. So by [Stein] 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

 (12) ∥(∇′)2uhi∥L∞(R2)+∥∇′nhi∥L∞(R2)≤1h,

while at the same time there is a set with

 (13) limsuph→0|Mh|h2=0

such that

 (14) uhi=ui and nhi=ni on S∖Mh.

Let and define the recovery sequence

 yh(x′,x3)=2∑i=1(uhi(x′)+hx3nhi(x′))ηhi(x′)+h2x232d(x′).

Since by (13), Lemma 3.3 shows that there is a constant depending only on and , such that choosing

 ρh=√2|Mh|δ

we have

 (15) Bρh(x0)∩S∖Mh≠∅ for all x0∈S.

By (13) we have

 (16) limsuph→0ρhh=0.

By (15) and (12), the fact that (hence ) almost everywhere on , implies that

 (17) 1h∥uhi−ui∥L∞+∥∇′uhi−∇′ui∥L∞+∥nhi−ni∥L∞≤Cρhh.

By (16) the right-hand side converges to zero. Hence (17) implies that and are bounded in since , .
Now we compute

 ∇′yh(x′,x3) =2∑i=1(∇′uhi(x′)+hx3∇′nhi(x′))ηhi(x′) +2∑i=1(uhi(x′)+hx3nhi(x′))∇′ηhi(x′)+h2x232∇′d(x′).

Recalling that and , we see that on

 (18) |∇′yh|≤2∑i=1(|∇′uhi|+h|∇′nhi|)+|uh1−uh2||∇′ηh|+h|nh1−nh2||∇′ηh|+h2|∇′d|≤C(1+1r|uh1−uh2|+hr|nh1−nh2|).

We have used the bound (11) as well as the uniform boundedness of and the estimate (12) for . Similarly, we can estimate

 (19) 1h|∂3yh|≤C(|nh1−nh2|+h).

Since the sequence remains uniformly bounded and recalling (7), we deduce from (18) and (19) that

 (20) |∇hyh|≤C(1+1r|uh1−uh2|) % on Σr×I.

We claim that

 (21) |uh1−uh2|≤Cr on Σr.

In fact, let . By definition of there exists such that . Since on and since is Lipschitz, we have

 |u1(x′)−u(ζ)|=|u1(x′)−u1(ζ)|≤C|x′−ζ|≤Cr.

Since and since is Lipschitz, we have

 |u1(x′)−u2(x′)|≤|u1(x′)−u(ζ)|+|u(x′)−u(ζ)|≤Cr.

After repeating the argument with and , we conclude that

 |u1−u2|≤Cr on Σr.

Together with (17) and (7), estimate (21) follows.

By (21) and (20) we see that is uniformly bounded in . Since is locally bounded, this implies

 limsuph→0∥W(∇hyh)∥L∞(Σr×I)<∞.

Therefore, since due to the regularity of ,

 1h2Eh(yh,Σr×I) =εh2∫Σr×IW(∇hyh)dx≤Cεrh2.

The right-hand side converges to zero due to (8).

On we have

 ∇hyh=(∇′uhi | nhi)+hx3(∇′nhi | d)+h2x232(∇′d | 0).

The map clearly takes values in . Define

 ahi=x3RT(∇′nhi | d)+hx232RT(∇′d | 0).

Estimates (12) and (17) ensure that

 (22) limsuph→0h∥ahi∥L∞(Si)<∞.

On we have , hence almost everywhere on this set. Hence

 (23) limsuph→0∥ahi∥L∞(Si∖Mh)<∞.

By the frame indifference of we have, almost everywhere on ,

 (24) W(∇hyh)=W