Processing math: 46%
Research article

Variational modeling of paperboard delamination under bending

  • Received: 11 February 2022 Revised: 15 May 2022 Accepted: 03 June 2022 Published: 15 June 2022
  • We develop and analyze a variational model for laminated paperboard. The model consists of a number of elastic sheets of a given thickness, which – at the expense of an energy per unit area – may delaminate. By providing an explicit construction for possible admissible deformations subject to boundary conditions that introduce a single bend, we discover a rich variety of energetic regimes. The regimes correspond to the experimentally observed: initial purely elastic response for small bending angle and the formation of a localized inelastic, delaminated hinge once the angle reaches a critical value. Our scaling upper bound then suggests the occurrence of several additional regimes as the angle increases. The upper bounds for the energy are partially matched by scaling lower bounds.

    Citation: Sergio Conti, Patrick Dondl, Julia Orlik. Variational modeling of paperboard delamination under bending[J]. Mathematics in Engineering, 2023, 5(2): 1-28. doi: 10.3934/mine.2023039

    Related Papers:

    [1] Francesco Maddalena, Danilo Percivale, Franco Tomarelli . Signorini problem as a variational limit of obstacle problems in nonlinear elasticity. Mathematics in Engineering, 2024, 6(2): 261-304. doi: 10.3934/mine.2024012
    [2] L. De Luca, M. Ponsiglione . Variational models in elasticity. Mathematics in Engineering, 2021, 3(2): 1-4. doi: 10.3934/mine.2021015
    [3] Dario Mazzoleni, Benedetta Pellacci . Calculus of variations and nonlinear analysis: advances and applications. Mathematics in Engineering, 2023, 5(3): 1-4. doi: 10.3934/mine.2023059
    [4] Francesco Freddi, Filippo Riva . Potential-based versus non potential-based cohesive models accounting for loading and unloading with application to sliding elastic laminates. Mathematics in Engineering, 2025, 7(3): 406-438. doi: 10.3934/mine.2025017
    [5] Federico Cluni, Vittorio Gusella, Dimitri Mugnai, Edoardo Proietti Lippi, Patrizia Pucci . A mixed operator approach to peridynamics. Mathematics in Engineering, 2023, 5(5): 1-22. doi: 10.3934/mine.2023082
    [6] Raimondo Penta, Ariel Ramírez-Torres, José Merodio, Reinaldo Rodríguez-Ramos . Effective governing equations for heterogenous porous media subject to inhomogeneous body forces. Mathematics in Engineering, 2021, 3(4): 1-17. doi: 10.3934/mine.2021033
    [7] Piermarco Cannarsa, Rossana Capuani, Pierre Cardaliaguet . C1;1-smoothness of constrained solutions in the calculus of variations with application to mean field games. Mathematics in Engineering, 2019, 1(1): 174-203. doi: 10.3934/Mine.2018.1.174
    [8] Manuel Friedrich . Griffith energies as small strain limit of nonlinear models for nonsimple brittle materials. Mathematics in Engineering, 2020, 2(1): 75-100. doi: 10.3934/mine.2020005
    [9] Luca Placidi, Emilio Barchiesi, Francesco dell'Isola, Valerii Maksimov, Anil Misra, Nasrin Rezaei, Angelo Scrofani, Dmitry Timofeev . On a hemi-variational formulation for a 2D elasto-plastic-damage strain gradient solid with granular microstructure. Mathematics in Engineering, 2023, 5(1): 1-24. doi: 10.3934/mine.2023021
    [10] Michiaki Onodera . Linear stability analysis of overdetermined problems with non-constant data. Mathematics in Engineering, 2023, 5(3): 1-18. doi: 10.3934/mine.2023048
  • We develop and analyze a variational model for laminated paperboard. The model consists of a number of elastic sheets of a given thickness, which – at the expense of an energy per unit area – may delaminate. By providing an explicit construction for possible admissible deformations subject to boundary conditions that introduce a single bend, we discover a rich variety of energetic regimes. The regimes correspond to the experimentally observed: initial purely elastic response for small bending angle and the formation of a localized inelastic, delaminated hinge once the angle reaches a critical value. Our scaling upper bound then suggests the occurrence of several additional regimes as the angle increases. The upper bounds for the energy are partially matched by scaling lower bounds.



    Paperboard is an important engineering material, widely used for packaging, e.g., in the food industry. Due to its sustainability, paper-based materials have more recently gained in interest also for other applications [22]. Paperboard is essentially a comparatively thick material made of processed wood pulp. The present article is concerned with laminated paperboard, which consists of multiple layers of paper bonded by an adhesive.

    The deformation of paperboard is a complex process taking place on multiple scales. Deep drawing of paperboard to achieve a desired geometry – similar to deep drawing of metals – is an important technique that is the subject of current research [18]. For an overview of current modeling approaches, from the individual fiber level to the laminate structures in laminated paperboard, see [22].

    Of particular interest is the formation of individual hinges in laminated paperboard when it undergoes bending. In [4,5,11,18,21], the formation of such hinges has been closely examined. As one can see in Figure 1, even in a simple experiment, such a hinge in laminated paperboard is characterized by a localized delamination of the individual sheets together with a localization of the deformation. A close up of the formed hinge is shown in Figure 2 (left). One can clearly see that the delaminated layers on the inside of the bend buckle, and thus release compressive stress. For experiments in a more controlled environment with an added crease to determine the exact location of the hinge, see, e.g., [4,Figure 3]. In particular in [4,5], the experimental observations are complemented by a finite element model, where the delamination is treated using cohesive zones. A variational model, together with a rigorous mathematical treatment, resulting in a rich phase diagram of different energetic scaling regimes for laminated paperboard undergoing a simple bend, is provided in this article.

    Figure 1.  Bending of laminated paperboard (BRAMANTE Buchbinderhartpappe 2mm) with increasing angle. a) unbent, b) purely elastic deformation, c, d) increasing delamination (red arrows), e), f) bending concentrates on delaminated hinge (red arrow). Images by D. Valainis.
    Figure 2.  Left: closeup of the hinge formed in the experiment shown in Figure 1f). One can clearly see the concentration of the bending angle to the hinge as well as the delaminated layers of paperboard. Image by D. Valainis. Right: hinge construction with energetically optimal scaling.
    Figure 3.  Piecewise affine construction. The angle formed by the segments with the x axis is α on the left, β in the first (downward) part of the central fold, then β, and finally α on the right.

    In our variational approach, the individual sheets are modeled using a nonlinearly elastic energy. The entire paperboard, but also of course the individual sheets of paper it consists of, are treated as slender elastic objects [2,3]. More specifically, the rigorous methods used to derive plate energies from nonlinear bulk elasticity [10,16,17] are used here to show lower scaling bounds for the deformed paperboard.

    Delamination is treated as Griffith-type fracture, with a fixed energetic cost per delaminated unit area [7,15], but the fracture surface is restricted to the interfaces of the individual sheets of paper in the laminated paperboard. We refer to [1] for an overview of the treatment of spaces of bounded variation which naturally occur in the mathematical treatment of fracture problems.

    The competition between these two energies, a bulk elastic energy and a surface fracture energy, leads to the emergence of different scaling regimes, depending on parameters. In this sense, our article is inspired by the seminal work of Kohn and Müller [14]. More closely related is the treatment of the blistering of thin films on a rigid substrate under compression [6,8,12], where the formation of self-similar, branching channels was observed. The scaling regimes observed in the present work should be compared to the energy scaling for the folding of single sheets of paper [9,19,23].

    Our main results can be summarized as follows. In our model, laminated paperboard exhibits different scaling regimes when subjected to bending boundary conditions, depending on the bending angle α, the paperboard thickness h and length L, number of sheets N and Griffith energy coefficient γ. For simplicity we focus here on the regime described by (3.74) below, which requires h>γN3, h5N5>γL4, and h5<γL4N3. Corresponding results for the other cases are discussed in Section 3, see in particular Remark 3.8 and Remark 3.9.

    First, for very small bending angle, no delamination occurs and a usual, energy minimizing single arch is formed. This purely elastic deformation, which corresponds to the bending of a thin plate, has an energy of order α2h3L.

    Then, sharply localized delamination occurs for all layers at a length scale α1/3h4/3γ1/3N. The energy scales as α1/3γ2/3h4/3.

    For even larger bending angle, the delamination length increases to αh3/2γ1/2N3/2 and the energy scales as αγ1/2h3/2N1/2.

    Finally, for larger bending angles, if the cost for delamination is small the entire paperboard may delaminate, which results in the standard bending energy for each of the N sheets individually, yielding a total energy of α2h3LN2.

    The occurrence of each individual regime depends on the magnitude of the model parameters. The second and third regime above correspond to the localized hinges observed in experiment. For most regimes, we provide rigorous lower bounds contingent on some assumptions on the delaminated sets. Explicit constructions realizing all energetic regimes by ensuring that individual delaminated sheets' midplanes deform isometrically are provided. These constructions introduce an additional fold on the inside of the hinges, which allows the individual delaminated sheets to deform isometrically, as seen in the buckling visible in Figure 2, where the experiment is shown on the left, and a construction with energetically optimal scaling is shown on the right.

    The remainder of this article is organized as follows. In Section 2 we introduce our mathematical model. Scaling upper and lower bounds for the energy are proved in Sections 3 and 4, respectively. We close with a brief discussion in Section 5.

    We consider a paperboard sample of thickness h, consistings of N layers, so that at most N1 delamination surfaces are possible. We let 2L be the length of the sample in the direction of bending, and consider a section of unit length in the third direction, which will be irrelevant for our results. We thus obtain a reference configuration Ωh:=(L,L)×(0,h)×(0,1).

    The set of admissible deformations consists of maps U:ΩhR3 which jump only on (a subset of) N1 prescribed planes which correspond to the delamination surfaces between layers, in the sense that the set JU of jump points of U obeys (up to null sets)

    JUΩh{x:x2hNZ}=(L,L)×{hN,2hN,,(N1)hN}×(0,1). (2.1)

    Outside its jump set the function is assumed to have a weak gradient which is square integrable. Mathematically, this means that U belongs to SBV2N(Ωh;R3), which we define as the space of functions in SBV(Ωh;R3) such that (2.1) holds and UL2(Ωh;R3×3). We recall that SBV(Ωh;R3) is the set of special functions of bounded variation, i.e., the set of integrable functions U:ΩhR3 such that the distributional gradient DU is a measure of the form DU=UL3+[U]νH2JU, with UL1(Ωh;R3×3), JU the H2-rectifiable jump set of U, [U]:JUR3 its jump, and ν:JUS2 the normal to JU. We refer to [1] for details of this definition and the relevant properties of the space.

    In order to introduce a boundary condition for the hinge, we prescribe the deformation gradient on both ends of the domain, i.e.,

    DU(x)=ˆRαforx1<L2,DU(x)=ˆRαforx1>L2, (2.2)

    where ˆRαSO(3) is a rotation of angle α around x3,

    ˆRα=(cosαsinα0sinαcosα0001). (2.3)

    The energy of a deformation USBV2N(Ωh;R3) consists of the sum of an elastic energy, depending on the absolutely continuous part of the deformation gradient U, and a delamination energy, which is proportional to the total area of the delaminated set JU. Specifically,

    E3Dh[U]:=ΩhW3D(U)dx+γH2(JU),forUSBV2N(Ωh;R3). (2.4)

    Here W3D:R3×3[0,) is an elastic energy density which obeys for some c>0

    1cdist2(ξ,SO(3))W3D(ξ)cdist2(ξ,SO(3))forallξR3×3. (2.5)

    The parameter γ>0 is the delamination energy per unit area, which in the present setting has the dimensions of a length, and H2 denotes the Hausdorff measure. In particular, H2(JU) is the area of the delaminated set.

    Remark 2.1. A number of comments are in order.

    1) If one assumes that the material is inhomogeneous, with moderate changes of the elastic energy W3D from layer to layer, the analysis is unchanged, up to a modification in the constant c in (2.5).

    2) The aim of this paper is to derive bounds from geometrically admissible delaminated states, using the Griffith-type delamination energy γH2(JU). More refined fracture models, for example including a cohesive zone, are beyond the scope of this paper.

    3) Although the manufacturing process of paperboard makes the material highly orthotropic, with different elastic properties in the lamination direction compared to in-plane, the key ingredient in the chosen model is the competition between the delamination energy and the effective bending stiffness. The introduced Griffith parameter γ can be seen as a delamination energy normalized by the bending stiffness, which in turn depends largely on the in-plane uniaxial compressive modulus.

    In this paper we assume that no structure arises in the x3-direction, in the sense that the deformation U takes the form

    U(x1,x2,x3)=u(x1,x2)+x3e3 (2.6)

    for some uSBV2N(ωh;R2), with ωh:=(L,L)×(0,h). The latter set is defined as the set of SBV functions such that uL2(ωh;R2×2) and

    Ju{xωh:x2hNZ}=(L,L)×{hN,2hN,,(N1)hN}, (2.7)

    and one easily sees that USBV2N(Ωh;R3) if and only if uSBV2N(ωh;R2).

    The energy then reduces to

    Eh[u]:=E3Dh[U]=(L,L)×(0,h)W2D(u)dx+γH1(Ju),foruSBV2N(ωh;R2), (2.8)

    where W2D is defined by W2D(ξ):=W3D(ξ+e3e3) for ξR2×2, and H1 denots the one-dimensional Hausdorff measure, which measures length. Here and below we identify ξR2×2 with the matrix ˆξR3×3 characterized by ˆξij=ξij for i,j=1,2, zero otherwise, and correspondingly for vectors. From (2.5) one easily obtains that the two-dimensional reduced energy obeys

    1cdist2(ξ,SO(2))W2D(ξ)cdist2(ξ,SO(2))forallξR2×2. (2.9)

    The boundary condition (2.2) in turn is equivalent to

    Du(x)=Rαforx1<L2,Du(x)=Rαforx1>L2, (2.10)

    where we denote by RαSO(2) the two-dimensional rotation of angle α,

    Rα:=(cosαsinαsinαcosα). (2.11)

    The objective is to find critical bending angle and energy bounds for the interplay of bending and delamination in terms of h, L, γ.

    We first consider the case that no delamination occurs. Then it is natural to expect that the deformation has constant curvature. We use here the classical construction for plate theories, with some simplifications which do not alter the scaling of the energy (for example, in (3.4) we simply use x2 as a factor, and not x2h2). We refer to [10] for a more general mathematical treatment.

    Lemma 3.1. For all h,L>0 and α[0,π2] there is a map uW1,(ωh;R2) which obeys the boundary condition (2.10) and such that

    Eh[u]cα2h3L. (3.1)

    The map u is injective, andW1,(ωh;R2)SBV2N(ωh;R2) for all N.

    Proof. We first fix an arc of circle for the central part,

    f(x1):=L2α(sin(2αx1/L)cos(2αx1/L)), (3.2)

    which obeys |f(x1)|=1 for all x1, and then extend it piecewise affine in the two boundary regions, setting

    v(x1):={f(L2)+(x1+L2)f(L2),ifL<x1<L2,f(x1),ifL2x1L2,f(L2)+(x1L2)f(L2),ifL2<x1<L. (3.3)

    One easily verifies that vW2,((L,L);R2), with |v|=1 and |v"|2α/L almost everywhere. We then define u:ωhR2 by

    u(x1,x2):=v(x1)+x2(v)(x1), (3.4)

    where (a1,a2):=(a2,a1) denotes counterclockwise rotation by 90 degrees. One easily verifies that uW1,(ωh;R2), Du(x1,x2)=Rα for ±x1(L2,L), and, for x1(L2,L2)

    Du(x)=v(x1)e1+(v)(x1)e2+x2(v")(x1)e1 (3.5)

    so that dist(Du(x),SO(2))|x2||v"(x1)|. Recalling the upper bound in (2.9) we obtain the desired bound

    Eh[u]LLh0c|x2v"(x1)|2dx2dx1=Lh0c(2αx2L)2dx2cα2h3L. (3.6)

    Injectivity can be easily verified from the definition of u.

    We start by constructing a continuous piecewise affine map ucpa that illustrates the basic structure of the deformation. This construction is not admissible for the functional considered here, but represents a deformation for the limiting case N. Before starting we introduce the notation

    Rφ:=(cosφsinφsinφcosφ)forφR. (3.7)

    We fix α(0,π2), h,L>0. We seek a map ucpa that obeys the boundary conditions (2.10), and that is symmetric with respect to the x2-axis, in the sense that

    ucpa(x1,x2)=(ucpa1(x1,x2)ucpa2(x1,x2)). (3.8)

    Therefore we can focus on {x10}, provided that we impose the condition ucpa1(0,x2)=0. We fix a parameter ζ(0,L/(2h)], and assume that

    Ducpa=Rαforζ(hx2)<x1<L,0<x2<h, (3.9)

    where Rα was defined in (2.11) (see Figure 3 for a sketch of the geometry). This ensures that the boundary condition (2.10) is fulfilled. In the region 0<x1<ζ(hx2) the material is allowed to shear and to open across the possible delamination lines, and of course to rotate by some angle β to be determined. However, to ensure global injectivity in the delaminated construction following in section 3.3, volumetric compression must be prevented. In other words we assume that the deformation gradient takes the form

    F:=Rβ(1a10a2) (3.10)

    for some a1R (representing shear) and a2[1,) (representing opening across the delamination lines). The symmetry condition (3.8) requires that ucpa1(0,x2)=0 for all x2(0,h). Since we already fixed ucpa1(0,h)=0, it suffices to require 1ucpa(0,x2)=0 for all x2(0,h), which is equivalent to F12=0 and therefore to a:=(a1,a2)=|a|(sinβ,cosβ). Denoting d:=|a| we obtain

    F=Rβ(1dsinβ0dcosβ)=(cosβ0sinβd) (3.11)

    with the condition dcosβ1. The construction is concluded if we impose continuity across the line {x1=ζ(hx2)}, which is equivalent to

    0=(RαF)(ζ1)=(ζcosαsinαζcosβζsinαcosαζsinβ+d). (3.12)

    In turn, this can be rewritten as the two conditions

    ζ=sinαcosαcosβ (3.13)

    and

    d=cosα+ζ(sinβ+sinα)=1cosαcosβ+sinαsinβcosαcosβ (3.14)

    which define ζ and d in terms of the angle β. It remains to verify the conditions a2=dcosβ1 and 0ζL/(2h), which limit the admissible choices of the angle β. We observe that this condition is needed to ensure injectivity of the final construction. We assume that α,β(0,π2); by ζ>0 we necessarily have β>α.

    The condition dcosβ1 can be rewritten as

    fα(β):=1cosαcosβ+sinαsinβcosαcosβcosβ1. (3.15)

    A short computation shows that

    fα(β)=sinα+β2sinβα2cosβ (3.16)

    and

    fα(β)=(cosβcosα)sinβcosβsinα2sin2αβ2<0 (3.17)

    so that for any α(0,π2) the function fα:(α,12π]R is strictly decreasing with fα(π2)=0 and fα(α+)=. Therefore we can define βeq(α)(α,π2) as the unique solution to fα(β)=1. By the implicit function theorem, βeqC1((0,π2)), and from

    fα(β)α=sinβcosβ1cos(αβ)>0 (3.18)

    we obtain βeq>0. The condition dcosβ1 is equivalent to β(α,βeq(α)]. We summarize and extend these results in the following statement.

    Lemma 3.2. (i) There is a continuous, increasing function βeq:(0,π2)(0,π2) such that α<βeq(α) for all α and

    1cosαcosβeq(α)+sinαsinβeq(α)cosαcosβeq(α)cosβeq(α)=1. (3.19)

    It obeys

    βeq(α)=41/3α1/3+o(α1/3)asα0. (3.20)

    (ii) For any α,β(0,π2), h,L>0, if

    α<ββeq(α)andζ:=sinαcosαcosβL2h (3.21)

    then the map ucpa:¯ωhR2 defined by

    ucpa(x):={(x1cosβ|x1|sinβ+dx2),if|x1|<ζ(hx2),(x1cosα+(x2h)sinαsgnx1|x1|sinα+(x2h)cosα+dh),ifζ(hx2)|x1|L, (3.22)

    is continuous, injective, obeys the boundary conditions and |1ucpa|=1 almost everywhere.

    Proof. (i) The function βeq is defined by fα(βeq(α))=1, with f as in (3.15). In order to verify (3.20), we write (3.19) as

    (1cosαcosβ+sinαsinβ)cosβ=cosαcosβ, (3.23)

    expand both sides to second order in α and rearrange terms, to obtain

    αsinβcosβ+12α2(cos2β+1)=(1cosβ)2+O(α3). (3.24)

    In particular, limα0βeq(α)=0. Expanding to leading order in β leads to αβ+α2=(β2/2)2+O(α3)+O(αβ3)+O(α2β2)+O(β6), which implies (3.20).

    (ii) We define ζ and d by (3.13) and (3.14), respectively. One verifies that ucpa as defined in (3.22) obeys the symmetry condition (3.8), ucpa1(0,x2)=0, and Ducpa=F for 0<x1<ζ(hx2), Ducpa=Rα for x1ζ(hx2), and that the function is continuous on the point (0,h). As is apparent from the construction above (see in particular (3.10)) this map is injective and transforms horizontal lines in a length-preserving way, in the sense that |1ucpa|=1. The rest follows from the computations above.

    The construction arises as a discretization and regularization of the continuous piecewise affine construction of Section 3.2. The final result of the construction, for two choices of β, is illustrated in Figure 4.

    Figure 4.  Construction of the multilayered folding construction. Top: β=βeq, bottom: β<βeq. This is a smoother version of the backbone structure from Figure 3. Also in this case the angle formed by the straight segments with the x axis is α, β, β and α (from left to right).

    We assume that in the inner region the material is partially delaminated, and use this to replace the non-isometric gradient F by a deformation which is isometric away from the discontinuity set. We then replace the sharp corners by regularized corners, in which each layer smoothly bends. The angle with the x1-axis is α at first, then gradually becomes β, then β, and finally α.

    We consider a construction with nN layers of paperboard which have been separated across n1 delamination surfaces. We assume n1; the case n=1 without delamination has already been treated in Section 3.1. To this end, we fix a subset {b1<b2<<bn1}(0,h)hNZ and assume that delamination occurs only on the surfaces {x2=b1,bn1}, in the sense that Ju(L,L)×{b1,,bn1}. For notational convenience we denote b0:=0 and bn:=h. We label by hj:=bj+1bj, 0j<n, the thickness of the j-th layer. The map we construct will be in C1([L,L]×[bj,bj+hj)) for each j<n.

    Step 1. The first step is a piecewise affine construction. We first construct the map on the set {x2=bj} using the continuous piecewise affine construction as background. Recalling (3.22) we define maps ˆf0,,ˆfn1:[L,L]R2 by

    ˆfj(x1):=ucpa(x1,bj)={(x1cosβ|x1|sinβ+dbj),if|x1|<lj,(x1cosα+(bjh)sinαsgnx1|x1|sinα+(bjh)cosα+dh),iflj|x1|L, (3.25)

    where lj:=ζ(hbj) and ζ is defined in (3.13). These maps are illustrated in Figure 3. The inner part of the construction, for |x1|lj, will be referred to as the "down slope".

    Step 2. In a next step we round the corners. This could be done by mollification, but it is important to (i) check the length, and keep the deformation isometric in the longitudinal direction, and (ii) be able to verify global injectivity of the two-dimensional deformation. Therefore we take a more explicit path and insert a circular arc of length arc>0 at each of the points of discontinuity of ˆfj, namely, at x1=±lj and at x1=0.

    In order to do this, we first give in (3.26) an explicit expression for the orientation φj(x1) of fj(x1) as a function of x1, and then in (3.27) construct fj, so that it is automatically parametrized by arc length. The position of the arcs is chosen so that, in the limit of small arc length arc0, we recover the piecewise affine function ˆfj, and that the boundary condition is fulfilled also for positive arc.

    The centers as well as radii are chosen such that the tangents of the arc match the respective left and right tangents of ˆfj, as illustrated in Figure 4. In particular, we can see that the radii are given by arcα+β and arcβ, respectively – independently of the layer index j. The requirement that the deformation is a rigid body motion for |x1|>L/2 results in the constraint that 2arc+max{lj}L/2, as neither the down-slope part nor the arcs can fulfill this condition.

    For each j, we first define the orientation function φjW1,(R) by

    φj(x1):={x1arcβ,if|x1|arc,βsgnx1,ifarc<|x1|<lj+arc,βsgnx1|x1|ljarcarc(α+β)sgnx1,iflj+arc|x1|lj+2arc,αsgnx1,if|x1|>lj+2arc, (3.26)

    and then define the arc-length preserving map of the j-th layer's mid-plane fj:[L,L]R2 by

    fj(x1):=ˆfj(0)+x10Rφj(x1)e1dx1. (3.27)

    We observe that fjW2,([L,L];R2), with |fj|=1 and |fj"|(α+β)/arc almost everywhere. By comparison with the derivative of ˆfj as defined in (3.25) we see that fj(0)=ˆfj(0), fj(x1+arc)=ˆfj(x1)=Rβe1 for x1(0,lj) and fj(x1+2arc)=ˆfj(x1)=Rαe1 for x1>lj. A short computation shows that if x1max{lj+2arc,lj+1+2arc} then

    fj+1(x1)fj(x1)=ˆfj+1(0)ˆfj(0)+(lj+1lj)(Rβe1Rαe1)=ˆfj+1(x1)ˆfj(x1)=Rαe2hj. (3.28)

    Further, we observe that

    ljl0=hζ=hsinαcosαcosβ. (3.29)

    Step 3. We are now ready to define the required mapping u:ωhR2. We use the same construction used in (3.4) in the proof of Lemma 3.1 with fj in place of v and set

    u(x1,x2):=fj(x1)+(x2bj)(fj)(x1),forx2[bj,bj+hj). (3.30)

    This map clearly belongs to SBV2N(ωh;R2).

    Assume now that l0+2arc12L. Then for x1L2 we have fj=Rαe1, and (3.28) holds. Therefore u is continuous for x1L2, with Du=Rα. The same holds (flipping some signs) on the other side, and therefore u fulfills the boundary condition (2.10).

    The energy is estimated by the same argument as in Section 3.1, see in particular (3.6). In particular, each of the 4n individual arcs in the construction has a change in angle of magnitude bounded by a universal constant times β, an arc-length no less than arc, and a thickness hj. In the rest of the domain the function fj is affine. Therefore

    Eh[u]n1j=12γ(lj+2arc)+n1j=0LLbj+1bjc|(x2bj)fj"(x1)|2dx2dx14nγ(hζ+arc)+cn1j=0β2h3jarc. (3.31)

    Step 4. We finally show that the map u defined in (3.30) is injective. This leads to an additional constraint.

    We first consider injectivity inside a single layer. For the affine parts and the arcs around ±(lj+arc) this follows by the same easy argument as in Lemma 3.1. For the central arc, with a different concavity, injectivity of the expression in (3.30) is equivalent to the fact that the layer thickness hj is not larger than the radius of the arc of circle described by fj, which is arc/β. Therefore injectivity is equivalent to

    hjarcβforj=0,,n1. (3.32)

    We next focus on the interaction between different layers. The boundary data automatically imply injectivity for |x1|12L; we can easily extend the construction to R×[0,h) and we see that it suffices to show that the set u(R×[bj,bj+hj)) does not intersect the curve (x1,fj+1(x1))=u(x1,bj+1). To ensure global injectivity we therefore need to show

    fj(x1)+λ(fj)(x1)fj+1(x1)foralljandx1,x1R,λ[0,hj). (3.33)

    Hence, it suffices to prove that

    |fj(s)fj+1(t)|hjforalljands,t. (3.34)

    We first consider the inner part of the construction. For |s|,|t|arc+lj+1, we have fj+1(t)=fj(t)+dhje2. We recall that |φj|β everywhere, so that (3.27) implies

    |e2(fj(s)fj(t))||e1(fj(s)fj(t))|tanβ. (3.35)

    Therefore, writing for brevity A:=fj(s)fj(t)R2,

    |fj+1(t)fj(s)|2=A21+(A2+dhj)2cos2βsin2βA22+(A22+2A2dhj+d2h2j)=1sin2βA22+2A2dhj+d2h2jd2h2jd2h2jsin2β=d2h2jcos2βh2j (3.36)

    where in the last step we used the condition dcosβ1 which follows from ββeq(α). For an illustration see Figure 5b.

    Figure 5.  Illustration of injectivity estimates. On the left, the location of subfigures a), b), and c) in the construction is shown. a) relevant trigonometric quantities for the outer part's construction, b) for the inner part, c) in area marked in red, both constructions are valid.

    For the outer part of the construction now let |s|,|t|arc. We argue similarly as for the inner part and for notational simplicity only consider t,s<0, the other side being a symmetric analog. Fix δ:=tan1(ζ) and orthogonal unit vectors ˉe1:=Rαδe1, ˉe2:=Rαδe2=(sin(δα)cos(δα)). From (3.28) and lj=lj+1+ζhj we have, for t+ζhjL/2, fj+1(t+ζhj)fj(t)=hjRα(e2+ζe1)=hj1+ζ2ˉe2. As fj+1(t+ζhj)=fj(t) for t+ζhj<arc, we conclude that in this range fj+1(t+ζhj)=fj(t)+hj1+ζ2ˉe2.

    We next show that

    |φj(x1)α+δ|δforallx1[L,0]. (3.37)

    To see this, note that the angle between the up-slope and ˉe1 is exactly equal to δ, where cosδ=hj1+ζ2hj. The angle ˉδ between the down-slope and ˉe1 satisfies cosˉδ=(dcosβ)hj1+ζ2hjcosδ, since ββeq and therefore dcosβ1. For an illustration, see Figure 5a. Alternatively, (3.37) can be proven algebraically. First, φj[β,α] for x10, so that it suffices to prove that p:=α+β2δ. We set m:=βα2 and express α and β in terms of m and p. As tanδ=ζ, by monotonicity of tan we need to check

    tanpζ=sinαcosαcosβ=12(1tanm1tanp). (3.38)

    On the other hand, by (3.16) the assumption ββeq (in the form fα(β)1) is the same as

    1tanm1+sin2psinpcosp=2tanp+1tanp. (3.39)

    As these two conditions are easily seen to be equivalent, (3.37) holds.

    Thus we have |ˉe2(fj(s)fj(t))||ˉe1(fj(s)fj(t))|tanδ for s,t0. Setting ˉA:=(ˉe1(fj(s)fj(t))ˉe2(fj(s)fj(t))), we calculate

    \begin{equation} \begin{split} |f_{j+1}(t+h_j\zeta)-f_j(s)|^2 & = \bar{A}_1^2 + (\bar{A}_2+\sqrt{1+\zeta^2}h_j)^2 \\ &\ge \left(\frac{1}{\tan^2\delta} +1\right) \bar{A}_2^2 + 2\bar{A}_2\sqrt{1+\zeta^2}h_j + (1+\zeta^2)h_j^2 \\ & = \frac{1+\zeta^2}{\zeta^2}\bar{A}_2^2 + 2\bar{A}_2\sqrt{1+\zeta^2}h_j + (1+\zeta^2)h_j^2 \\ & \ge h_j^2. \end{split} \end{equation} (3.40)

    For the mixed case where s , t are in different parts of the we note that the minimal distance between two consecutive layers is achieved (again, on the left side of the folding construction) between diagonal corners of the red quadrilateral in Figure 5c, as, within this area, both the inner and the outer estimate are valid. It is clear that the respective distances are bounded by \min\{dh, \sqrt{1+\zeta^2}h\} .

    Lemma 3.3. Fix \alpha\in (0, \frac\pi4] , h > 0 , L \ge h , N\in \mathbb N , N\ge 1 . For any n\in \mathbb N with 1\le n\le N , any \beta\in (\alpha, \beta_{{\mathrm{eq}}}(\alpha)] , and any \ell_{\mathrm{arc}} > 0 , if

    \begin{equation} \frac{2\beta h}{n} \le \ell_{\mathrm{arc}}\le \frac 18 L \hskip5mm \mathit{and}\hskip5mm \ell: = \frac{\sin\alpha}{\cos\alpha-\cos\beta}h \le\frac{1}{4}L \end{equation} (3.41)

    then there exists a map u\in SBV^2_N(\omega_h; \mathbb R^2) , which obeys (2.10), is injective, and such that

    \begin{equation} {\cal E}_h[u]\le c\left(\gamma (\ell+ \ell_{\mathrm{arc}})n + \frac{\beta^2 h^3}{ \ell_{\mathrm{arc}} n^2}\right) = c\left(\gamma n\left(\frac{ h\sin \alpha}{\cos\alpha - \cos\beta}+ \ell_{\mathrm{arc}}\right) + \frac{\beta^2 h^3}{ \ell_{\mathrm{arc}} n^2 }\right) \end{equation} (3.42)

    with a constant c > 0 only depending on the elastic energy density W_{\mathrm{2D}} .

    Proof. For n = 1 there is no delamination, the assertion follows from Lemma 3.1, \alpha\le\beta , and \ell_{\mathrm{arc}}\le L . Therefore we can assume n\ge 2 in the following.

    We first choose a subset \{b_1, \dots, b_{n-1}\}\subseteq \{h\frac 1N, \dots, h\frac{N-1}N\} such that, setting b_0 = 0 and b_{n} = h , we have h_j = b_{j+1}-b_j\le \frac{2h}{n} for all j = 0, \dots, n-1 . We immediately notice that the map u given in (3.30) is a member of SBV^2_N(\omega_h; \mathbb R^2) , with its jump set contained in the set (-2 \ell_{\mathrm{arc}} - \ell, 2 \ell_{\mathrm{arc}} + \ell)\times\{b_1, \dots, b_{n-1}\} . The energy estimate follows then from (3.31).

    The conditions in (3.41) imply in particular that 2 \ell_{\mathrm{arc}}+\ell\le L/2 , and therefore that the boundary condition (2.10) is fulfilled. The condition (3.32), required for injectivity around the central arc, follows from the first inequality in (3.41) and the condition h_j\le 2h/n .

    Remark 3.4. In reality we expect that buckling only appears in the lower layers in Figure 4, while the upper layer may still be under small tension. A more complex construction taking this into account would, however, only yield a constant factor in the energy and is therefore not needed in order to derive the scaling bounds.

    Theorem 3.5. There are C > 0 and \eta\in(0, 1] such that for 0 < h\le \eta L , \alpha\in (0, \pi/2] there is u\in SBV^2_N(\Omega; \mathbb R^2) which obeys (2.10), u\in W^{1, \infty}(\Omega\setminus\overline{J_u}; \mathbb R^2) with \mathcal H^1(\overline{J_u}\setminus {J_u}) = 0 , is injective, and obeys

    \begin{equation} {\cal E}_h[u]\le C \min\{ \frac{\alpha^2h^3}{L}, \alpha^{1/3} \gamma^{2/3}h^{4/3} +\frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}} + \frac{\alpha^2 h^3}{L N^2 } { + \frac{\alpha h^4}{L^2 N^2 }}\}. \end{equation} (3.43)

    Proof. The result in Theorem 3.5 will follow from Propositions 3.6 and 3.7. The first proposition considers the case where the angle \alpha is large, the second treats the case of small \alpha . As \frac{ h^3}{L N^2 }\le \frac{h^4}{L^2 N^2 } , in the case of large \alpha the last term does not appear.

    Proposition 3.6. There are C > 0 and \eta\in(0, 1) such that for 0 < h\le \eta L , \alpha\in (0, \pi/2] there is u\in SBV^2_N(\Omega; \mathbb R^2) which obeys (2.10), is W^{1, \infty} in \Omega\setminus\overline{J_u} , injective, and obeys

    \begin{equation} {\cal E}_h[u]\le C \min\{ \frac{h^3}{L}, \gamma^{2/3}h^{4/3} +\frac{ \gamma^{1/2} h^{3/2}}{N^{1/2}} + \frac{ h^3}{L N^2 }\}. \end{equation} (3.44)

    Proof. We prove the proposition in 3 steps.

    Step 1. If \alpha > \frac\pi4 we can combine two folds with L': = L/4 and \alpha': = \alpha/2 . Therefore we assume in the remainder of the proof that \alpha\in(0, \frac\pi4] .

    Step 2. By Lemma 3.1 there is an injective map u\in W^{1, \infty}(\omega_h; \mathbb R^2) with {\cal E}_h[u]\le C h^3/ L . This proves the first bound.

    Step 3. We intend to use Lemma 3.3 with \beta: = \beta_{{\mathrm{eq}}}(\alpha) . We first check the second condition in (3.41). The function

    \begin{equation} \alpha\mapsto \frac{\sin\alpha}{\cos\alpha-\cos \beta_{{\mathrm{eq}}}(\alpha)} \end{equation} (3.45)

    is continuous on (0, \frac\pi4] , and as \beta_{{\mathrm{eq}}} = \mathcal O(\alpha^{1/3}) for small \alpha , it converges to 0 for \alpha\to0 . Therefore there is C > 0 such that

    \begin{equation} 0 < \frac{\sin\alpha}{\cos\alpha-\cos \beta_{{\mathrm{eq}}}(\alpha)}\le C \hskip5mm \;{\rm{ for\; all }}\; \alpha\in(0,\frac{\pi}{4}]. \end{equation} (3.46)

    By Lemma 3.3 we thus know that for any n\in\{1, \dots, N\} and \ell_{\mathrm{arc}} > 0 , if

    \begin{equation} \frac{2\beta h}{n} \le \ell_{\mathrm{arc}}\le \frac 18 L \hskip5mm \;{\rm{ and }}\;\hskip5mm C h \le\frac{1}{4}L \end{equation} (3.47)

    then

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le c\left( \gamma h n+\gamma \ell_{\mathrm{arc}} n+ \frac{ h^3}{ \ell_{\mathrm{arc}}\; n^2 }\right). \end{equation} (3.48)

    We assume that

    \begin{equation} \frac hL \le \min\left\{\frac{1}{4C},\frac1{8\pi}\right\}, \end{equation} (3.49)

    which is guaranteed if \eta is chosen appropriately. Noting that reducing \ell_{\mathrm{arc}} below \mathcal O(h) does not reduce the energy, due to the first term in (3.48), we can choose

    \begin{equation} \ell_{\mathrm{arc}}\in \left[\pi h, \frac L8\right]. \end{equation} (3.50)

    Since \beta\le \frac\pi2 , we are automatically assured that all requirements in (3.47) are satisfied. At the same time, if \ell_{\mathrm{arc}} is chosen in this range then we can ignore the first term in (3.48). Condition (3.49) ensures that the set of possible values of \ell_{\mathrm{arc}} is nonempty.

    We first optimize the choice of n . Set

    \begin{equation} n: = \min\left\{N,\left\lceil \frac{h}{\gamma^{1/3}\; \ell_{\mathrm{arc}}^{2/3}}\right\rceil\right\}\in\{1, \dots, N\}. \end{equation} (3.51)

    Inserting in (3.48) leads to

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le \begin{cases} c\left(\gamma \ell_{\mathrm{arc}} + \gamma^{2/3}h \ell_{\mathrm{arc}}^{1/3}\right) ,& \;{\rm{ if }}\; h\le N \gamma^{1/3} \ell_{\mathrm{arc}}^{2/3},\\ c \frac{ h^3}{ \ell_{\mathrm{arc}}\; N^2 }, & \;{\rm{ otherwise,}}\; \end{cases} \end{equation} (3.52)

    where the term \gamma \ell_{\mathrm{arc}} stems from the fact the n\ge 1 .

    It remains to choose the value of \ell_{\mathrm{arc}} . Since the first expression is strictly increasing, and the second one strictly decreasing, if the critical value is admissible then it is optimal. Precisely, if \ell_{\mathrm{arc}}^\; {\rm{crit}}\; : = \gamma^{-1/2} h^{3/2}N^{-3/2} \in[\pi h, L/8] then we obtain a bound of the form h^{3/2}\gamma^{1/2}/N^{1/2} . Otherwise, the relevant constraint is the one that prevents this optimal value from being chosen. Specifically,

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le \begin{cases} c\left(\gamma h + \gamma^{2/3}h^{4/3}\right) ,& \;{\rm{ if }}\; \gamma^{-1/2} h^{3/2}N^{-3/2} < \pi h,\\ c\frac{\gamma^{1/2}h^{3/2}}{N^{1/2}}, &\;{\rm{ if }}\; \gamma^{-1/2} h^{3/2}N^{-3/2}\in[\pi h, L/8],\\ c \frac{ h^3}{L N^2 }, & \;{\rm{ if }}\; L < 8 \gamma^{-1/2}h^{3/2}N^{-3/2}, \end{cases} \end{equation} (3.53)

    or, equivalently,

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le \begin{cases} c\left(\gamma h + \gamma^{2/3}h^{4/3}\right) ,& \;{\rm{ if }}\; h < \gamma N^3,\\ c\frac{\gamma^{1/2}h^{3/2}}{N^{1/2}}, &\;{\rm{ if }}\; \gamma N^3\le h\le \gamma^{1/3}L^{2/3}N,\\ c \frac{ h^3}{L N^2 }, & \;{\rm{ if }}\; \gamma^{1/3}L^{2/3}N < h. \end{cases} \end{equation} (3.54)

    We remark that the first and the last regime are disjoint, since \gamma^{1/3}L^{2/3}N < h and h\le L imply \gamma^{1/3}h^{2/3}N < h , which is equivalent to \gamma N^3 < h .

    The calculation above thus reveals energetic regimes as follows, which differ in the number of delaminated layers and the delamination length.

    Sharp fold partial delamination: The energy scaling \gamma h + \gamma^{2/3}h^{4/3} has the shortest possible delamination length \ell_{\mathrm{arc}} = \pi h and n\sim \lceil \gamma^{-1/3}h^{1/3}\rceil . It originates by balancing the two terms \gamma \ell_{\mathrm{arc}} n+ \frac{ h^3}{ \ell_{\mathrm{arc}} n^2 } after setting \ell_{\mathrm{arc}} = h .

    The term \gamma h , corresponding to n = 1 , can be dropped. To see this, we first note that it is only relevant if h\le\gamma . In this regime, however, it would be convenient not to delaminate at all, obtaining an energy h^3 L^{-1} . Indeed, as h\le L , h\le\gamma implies h^3 L^{-1} \le h^2 \le \gamma h .

    Localized full delamination: The energy scales as \frac{\gamma^{1/2}h^{3/2}}{N^{1/2}} . As n = N , each layer is delaminated, however only over a length \ell_{\mathrm{arc}} = \gamma^{-1/2}h^{3/2}N^{-3/2} . This originates by balancing the two terms \gamma \ell_{\mathrm{arc}} n+ \frac{ h^3}{ \ell_{\mathrm{arc}} n^2 } after setting n = N . All layers are delaminated, but only over a length \ell_{\mathrm{arc}} . The energy corresponds to N plates of thickness h/N bent over a length \ell_{\mathrm{arc}} , the value is determined so that this exactly balances the delamination energy.

    Total delamination: The energy scales as \frac{ h^3}{L N^2 } , there are N separate plates of thickness h/N , each bending over the entire available length L , the delamination energy is smaller.

    Since the powers of h are increasing, and the expression is continuous, one can also write the above regimes concisely as in (3.44).

    Proposition 3.7. There are C > 0 , \alpha_0\in(0, \pi/2] such that for 0 < h\le L , \alpha\in (0, \alpha_0] there is u\in SBV^2_N(\Omega; \mathbb R^2) which obeys (2.10), is W^{1, \infty} in \Omega\setminus\overline{J_u} , injective, and obeys

    \begin{equation} {\cal E}_h[u]\le C \min\left\{ \frac{\alpha^2 h^3}{L}, \alpha^{1/3} \gamma^{2/3}h^{4/3} +\frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}} + \frac{\alpha^2 h^3}{L N^2 } + \frac{\alpha h^4}{L^2 N^2 } \right\}. \end{equation} (3.55)

    Proof. As above, we prove the two bounds separately.

    Step 1. Again, we use Lemma 3.1 to obtain an injective map u\in W^{1, \infty}(\omega_h; \mathbb R^2) with {\cal E}_h[u]\le C \alpha^2h^3/ L .

    Step 2. The delaminated construction is obtained with Lemma 3.3 with a careful choice of the parameters. In this regime of small \alpha , any delaminated construction must also satisfy that \beta is small, as \beta \le \beta_{{\mathrm{eq}}}(\alpha)\simeq (4\alpha)^{1/3} for \alpha\to0 . A straightforward Taylor series expansion of the delamination length factor \zeta leads to

    \begin{equation} \zeta = \frac{\sin\alpha}{\cos\alpha-\cos\beta} = \frac{2\alpha + O(\alpha^3)} {\beta^2-\alpha^2+ O(\beta^4)+O(\alpha^4)}. \end{equation} (3.56)

    Since \alpha\le \beta , the O(\alpha^4) in the denominator can be ignored. It is however important to avoid that \beta^2-\alpha^2 in the denominator becomes too small. For this reason we shall restrict the choice of \beta by assuming \beta\ge 2\alpha , which implies \cos\alpha-\cos\beta\ge \cos\frac\beta2-\cos\beta = \frac38 \beta^2+O(\beta^4) . For \beta sufficiently small, this is larger than \beta^2/3 . Therefore there is \alpha_0 > 0 such that for all \alpha\in (0, \alpha_0] one has 2\alpha\le \beta_{{\mathrm{eq}}}(\alpha) and for all \beta\in [2\alpha, \beta_{{\mathrm{eq}}}(\alpha)]

    \begin{equation} \zeta = \frac{\sin\alpha}{\cos\alpha-\cos\beta}\le 3\frac{\alpha } {\beta^2}. \end{equation} (3.57)

    By Lemma 3.3 we know that for any n\in\{1, \dots, N\} and \ell_{\mathrm{arc}} > 0 , if

    \begin{equation} \frac{2\beta h}{n} \le \ell_{\mathrm{arc}} \le \frac 18 L \hskip5mm \;{\rm{ and }}\;\hskip5mm 3 \frac{\alpha}{\beta^2}h \le\frac{1}{4}L \end{equation} (3.58)

    (noting that h\le L guarantees that the set of admissible \beta\le \beta_{{\mathrm{eq}}} is not empty, after potentially decreasing \alpha_0 again) then

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le c\left( \frac{\alpha}{\beta^2}\gamma h n +\gamma \ell_{\mathrm{arc}} n+ \frac{\beta^2 h^3}{ \ell_{\mathrm{arc}}\; n^2 }\right). \end{equation} (3.59)

    We assume that

    \begin{equation} \frac{\alpha}{\beta^2}h\le \ell_{\mathrm{arc}}\le \frac18 L \end{equation} (3.60)

    so that the first term in the energy can be ignored. We next compare the first lower bound in (3.58) with the one in (3.60). Since we know that \beta\le \beta_{{\mathrm{eq}}}(\alpha)\le c \alpha^{1/3} , and n\ge 1 , we obtain for some c_* > 0

    \begin{equation} \frac{2\beta h}{n}\le c_*\frac{\alpha}{\beta^2}h. \end{equation} (3.61)

    Without loss of generality we can assume c_*\ge\frac32 . We then see that all conditions in (3.58) and in (3.60) are satisfied provided

    \begin{equation} c_*\frac{\alpha}{\beta^2}h\le \ell_{\mathrm{arc}}\le \frac18 L. \end{equation} (3.62)

    By the same argument we see that this condition is optimal (up to a factor). We choose

    \begin{equation} n: = \min\left\{N,\left\lceil \frac{\beta^{2/3}h}{\gamma^{1/3} \;\ell_{\mathrm{arc}}^{2/3}}\right\rceil\right\}\in\{1, \dots, N\}. \end{equation} (3.63)

    Inserting in (3.59) leads to

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le \begin{cases} c\left(\gamma \ell_{\mathrm{arc}} + \beta^{2/3} \gamma^{2/3} h \ell_{\mathrm{arc}}^{1/3}\right) ,& \;{\rm{ if }}\; \beta^{2/3}h\le \gamma^{1/3} \ell_{\mathrm{arc}}^{2/3} N,\\ c \frac{ \beta^2 h^3}{ \ell_{\mathrm{arc}} N^2 }, & \;{\rm{ otherwise,}}\; \end{cases} \end{equation} (3.64)

    noting again that the term \gamma \ell_{\mathrm{arc}} stems from rounding n up to the next integer, which is needed as we require n\ge1 . As in Step 3 of the proof of Proposition 3.6, the first expression is increasing in \ell_{\mathrm{arc}} , the second decreasing. Therefore the optimal value is the critical one, \ell_{\mathrm{arc}}^\; {\rm{crit}}\; : = \beta \gamma^{-1/2} h^{3/2}N^{-3/2} , if it is admissible in the sense of (3.62). If not, the optimal value is the admissible value closest to the critical one.

    Assume now that \ell_{\mathrm{arc}} = \ell_{\mathrm{arc}}^\; {\rm{crit}}\; is admissible. Then one obtains n = N , and \inf_u {\cal E}_h[u]\le c\frac{\beta \gamma^{1/2}h^{3/2}}{N^{1/2}} . Restating the requirements on \ell_{\mathrm{arc}} from (3.62), this is possible if \ell_{\mathrm{arc}}^\; {\rm{crit}}\; \ge c_*\alpha h/\beta^{2} and \ell_{\mathrm{arc}}^\; {\rm{crit}}\; \le L/8 , which is the same as

    \begin{equation} h \ge c_*^2\alpha^2\beta^{-6}\gamma N^3 \end{equation} (3.65)

    and

    \begin{equation} h \le \frac14 \beta^{-2/3}\gamma^{1/3}L^{2/3} N . \end{equation} (3.66)

    Assume now that one of these two conditions is violated (one can easily check that if c_*\frac{\alpha}{\beta^2}h\le \frac18 L , so that (3.62) gives a nonempty set of admissible values of \ell_{\mathrm{arc}} , then it is not possible for both (3.65) and (3.66) to be violated at the same time).

    ● If (3.65) is violated, then \ell_{\mathrm{arc}} = c_* \alpha \beta^{-2}h , n = \left\lceil \frac{\beta^{2/3}h}{\gamma^{1/3} \;\ell_{\mathrm{arc}}^{2/3}}\right\rceil , and \inf_u {\cal E}_h[u]\le c(\alpha\beta^{-2}\gamma h +\alpha^{1/3}\gamma^{2/3} h^{4/3}) .

    ● If (3.66) is violated, then \ell_{\mathrm{arc}} = L/8 , n = N , and \inf_u {\cal E}_h[u]\le c \frac{ \beta^2h^3}{L N^2 } .

    Summarizing, this leads to

    \begin{equation} E: = \inf\limits_u {\cal E}_h[u]\le \begin{cases} c\left(\frac{\alpha\gamma h}{\beta^{2} } + \alpha^{1/3}\gamma^{2/3}h^{4/3}\right) ,& \;{\rm{ if }}\; h < c_*^2\alpha^2\beta^{-6}\gamma N^3 ,\\ c\frac{\beta \gamma^{1/2}h^{3/2}}{N^{1/2}}, &\;{\rm{ if }}\; c_*^2\alpha^2\beta^{-6}\gamma N^3 \le h \le \frac14 \beta^{-2/3}\gamma^{1/3}L^{2/3} N,\\ c \frac{ \beta^2h^3}{L N^2 }, & \;{\rm{ if }}\; h > \frac14 \beta^{-2/3}\gamma^{1/3}L^{2/3} N, \end{cases} \end{equation} (3.67)

    or, equivalently,

    \begin{equation} E \le \begin{cases} c\left( \frac{\alpha\gamma h}{\beta^{2} } + \alpha^{1/3} \gamma^{2/3}h^{4/3}\right) ,& \;{\rm{ if }}\; \beta < \alpha^{1/3}\gamma^{1/6}h^{-1/6}N^{1/2},\\ c\frac{\beta \gamma^{1/2}h^{3/2}}{N^{1/2}}, &\;{\rm{ if }}\; \alpha^{1/3}\gamma^{1/6}h^{-1/6}N^{1/2}\le\beta\le \gamma^{1/2} h^{-3/2}L N^{3/2}, \\ c \frac{ \beta^2h^3}{L N^2 }, & \;{\rm{ if }}\; \beta > \gamma^{1/2}h^{-3/2}LN^{3/2}. \end{cases} \end{equation} (3.68)

    It remains to choose \beta\in[2\alpha, \beta_{{\mathrm{eq}}}(\alpha)] with c_*\alpha\beta^{-2}h\le \frac18 L (from (3.62)), i.e.,

    \begin{equation} \max\{2\alpha, (8c_*\alpha h/L)^{1/2}\}\le\beta\le \beta_{{\mathrm{eq}}}(\alpha). \end{equation} (3.69)

    These bounds scale as \alpha+\alpha^{1/2}(h/L)^{1/2} and \alpha^{1/3} , respectively, and in particular, possibly reducing \alpha_0 , we can ensure that the set is not empty for all values of the parameters.

    We observe that the expression in (3.68) is minimized at \beta_\; {\rm{crit}}\; : = \alpha^{1/3}\gamma^{1/6}h^{-1/6}N^{1/2} , and that – provided \alpha_0 is chosen sufficiently small – \alpha^{1/3}\le \beta_{{\mathrm{eq}}}(\alpha) for all \alpha\in(0, \alpha_0] . There are different regimes, both of which are characterized by the delamination length being at the minimum value still satisfying the injectivity conditions. In these two regimes the optimal scaling is attained by the choice \beta = \beta_{{\mathrm{eq}}}(\alpha) , which corresponds to the fact that the layers touch each other, as illustrated in Figure 4. The first one delaminates some of the layers, the second one all layers.

    Sharp fold partial delamination: If \beta_\; {\rm{crit}}\; > \alpha^{1/3} , which is the same as h < \gamma N^3 , then we take \beta = \alpha^{1/3} , and we are in the first regime in (3.68). This gives E\le c\alpha^{1/3}\gamma h + c\alpha^{1/3} \gamma^{2/3}h^{4/3} . The first term does not contribute. Indeed, if \alpha^{1/3} \gamma^{2/3}h^{4/3}\le \alpha^2 h^3/L then necessarily \gamma^{2/3}L\le \alpha^{5/3}h^{5/3} , which, since h\le L and \alpha\le 1 , implies \gamma\le h . Therefore we can drop the \alpha^{1/3}\gamma h term and obtain E\le c\min\{\alpha^{1/3} \gamma^{2/3}h^{4/3}, \alpha^2 h^3/L\} . In this regime \ell_{\mathrm{arc}} = \alpha^{1/3}h and n\simeq \gamma^{-1/3}h^{1/3}\le N , so that the energy bound behaves as \gamma \ell_{\mathrm{arc}} n .

    Sharp fold full delamination: If \max\{2\alpha, (8c_*\alpha h/L)^{1/2}\}\le \beta_\; {\rm{crit}}\; \le \alpha^{1/3} , then \beta = \beta_\; {\rm{crit}}\; is optimal and E\le c\alpha^{1/3}\gamma^{2/3}h^{4/3} . In this regime n = N and \ell_{\mathrm{arc}} = \alpha^{1/3}\gamma^{-1/3}h^{4/3} N^{-1} (indeed, the bound on E is of the form \gamma \ell_{\mathrm{arc}} N ).

    If \beta_\; {\rm{crit}}\; < \max\{2\alpha, (8c_*\alpha h/L)^{1/2}\} then we take \beta = \max\{2\alpha, (8c_*\alpha h/L)^{1/2}\} , and we are in the second or third regime in (3.68). This gives in principle four subcases, one of which cannot occur. The first of the three possible cases still delaminates only locally, the other two are delaminated over the entire length of the specimen. In any case, all layers are delaminated.

    Localized full delamination: Assume first that \alpha\ge 2c_* h/L , so that \beta = 2\alpha . If 2\alpha\le\gamma^{1/2}h^{-3/2}LN^{3/2} then we are in the second regime in (3.68) (with n = N and \ell_{\mathrm{arc}} = \ell_{\mathrm{arc}}^\; {\rm{crit}}\; = 2\alpha \gamma^{-1/2} h^{3/2}N^{-3/2} ) and E\le c\frac{\alpha \gamma^{1/2}h^{3/2}}{N^{1/2}} .

    Total delamination: If \alpha\ge 2c_* h/L and 2\alpha > \gamma^{1/2}h^{-3/2} L N^{3/2} , then we are in the third regime in (3.68) (with n = N and \ell_{\mathrm{arc}} = L/8 ) and E\le c\frac{ \alpha^2h^3}{L N^2 } .

    Small-angle total delamination: If instead \alpha < 2c_* h/L then \beta = (8c_*\alpha h/L)^{1/2} . A short computation shows that \beta_\; {\rm{crit}}\; \le (\alpha h/L)^{1/2} is equivalent to (\alpha h/L)^{1/2}\ge \gamma^{1/2}h^{-3/2}LN^{3/2} , so that we are necessarily in the third regime in (3.68) (with n = N and \ell_{\mathrm{arc}} = L/8 ) and E\le c \frac{\alpha h^4}{ L^{2}N^{2}} .

    In all cases we may conclude

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le c\min\left\{ \frac{\alpha^2h^3}{L}, \alpha^{1/3}\gamma^{2/3}h^{4/3} +\frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}} {+ \frac{\alpha h^4}{L^2 N^2 }} + \frac{\alpha^2 h^3}{L N^2 } \right\}. \end{equation} (3.70)

    Remark 3.8. A more precise summary of the cases above shows that for h\le \gamma N^3 , then only the "localized partial delamination" regime is relevant, and one obtains (for \alpha\in(0, \alpha_0] )

    \begin{equation} \inf\limits_u {\cal E}_h[u]\le c\min\left\{ \frac{\alpha^2h^3}{L}, \alpha^{1/3}\gamma^{2/3}h^{4/3} \right\}. \end{equation} (3.71)

    Recalling (3.54) one can see that this holds for all \alpha\in(0, \frac\pi2] . As the delamination length is proportional to \alpha^{1/3}h , it is always smaller than L . In turn, the number of delaminated layers n\simeq \gamma^{-1/3}h^{1/3} does not depend on the opening \alpha . Both the delamination length and the number of delaminated layers are discontinuous at yielding point \alpha\simeq \gamma^{2/5}L^{3/5}/h .

    Remark 3.9. The statement in Theorem 3.5 (analogously in Propositions 3.6 and 3.7), can also be written in the form

    \begin{equation} {\cal E}_h[u]\le C \min\left\{ \frac{\alpha^2h^3}{L}, \max\bigg\{\alpha^{1/3} \gamma^{2/3}h^{4/3} ,\frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}} , { \frac{\alpha h^4}{L^2 N^2 }} ,\frac{\alpha^2 h^3}{L N^2 } \bigg\}\right\}. \end{equation} (3.72)

    To see this, it suffices to observe that \max\{a, b\}\sim a+b for any a, b\ge0 .

    In order to gain some understanding in these many regimes we sort them in increasing values of \alpha (see also Figure 6 below). For definiteness let us assume that the first linear term is the relevant one, in the sense that \gamma^{1/2}h^{3/2}N^{-1/2}\ge h^4L^{-2}N^{-2} , so that the last term in the maximum can be ignored. The condition is equivalent to h^5\le \gamma L^4 N^3 . A simple computation shows that the maximum is then equal to

    \begin{equation} \begin{cases} \alpha^{1/3} \gamma^{2/3}h^{4/3}, &\;{\rm{ if }}\; \alpha < \gamma^{1/4}N^{3/4}/h^{1/4},\\ \frac{\alpha \gamma^{1/2}h^{3/2}}{N^{1/2}},&\;{\rm{ if }}\; \gamma^{1/4}N^{3/4}/h^{1/4}\le\alpha\le \gamma^{1/2}LN^{3/2} /h^{3/2},\\ \frac{\alpha^2 h^3}{L N^2 },&\;{\rm{ if }}\; \alpha > \gamma^{1/2} L N^{3/2} /h^{3/2}. \end{cases} \end{equation} (3.73)
    Figure 6.  Left: sketch of the upper bound of the energy from Theorem 3.5, depending on the bending angle \alpha , in the regime of (3.74). Right: sketch of the moment (i.e., the derivative of the energy) depending on the bending angle \alpha . In both figures, the colors indicate the respective regime (blue: non-delaminated, red: localized full delamination. yellow: localized full delamination redux, purple: total delamination). The plots are a representation of an energy with the respective scaling regimes. Multiplicative and additive constants have been applied to produce a graph with continuity properties in line with physical considerations. The parameters for the plots displayed here are h = 1 , L = 10 , N = 8 , \gamma = 10^{-6} .

    It remains to take the minimum between this expression and \alpha^2h^3/L . This is always larger than the third expression, and one easily computes the points where this intersects the other two. There are two cases. If \gamma^{2/5}L^{3/5}/h < \gamma^{1/4}N^{3/4}/h^{1/4} (which is the same as \gamma L^4 < h^5N^5 ) then we obtain

    \begin{equation} \frac1c {\cal E}_h[u]\le \begin{cases} \frac{\alpha^2h^3}{L}, & \;{\rm{ if }}\; \alpha < \gamma^{2/5}L^{3/5}/h,\\ \alpha^{1/3} \gamma^{2/3}h^{4/3}, &\;{\rm{ if }}\; \gamma^{2/5}L^{3/5}/h\le\alpha < \gamma^{1/4}N^{3/4}/h^{1/4},\\ \frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}},&\;{\rm{ if }}\; \gamma^{1/4}N^{3/4}/h^{1/4}\le\alpha < \gamma^{1/2} L N^{3/2} /h^{3/2},\\ \frac{\alpha^2 h^3}{L N^2 },&\;{\rm{ if }}\; \alpha\ge \gamma^{1/2}L N^{3/2} /h^{3/2}. \end{cases} \end{equation} (3.74)

    whereas for \gamma^{2/5}L^{3/5}/h\ge \gamma^{1/4}N^{3/4}/h^{1/4} the \alpha^{1/3} regime disappears and

    \begin{equation} \frac1c {\cal E}_h[u]\le \begin{cases} \frac{\alpha^2h^3}{L}, & \;{\rm{ if }}\; \alpha < \gamma^{1/2}LN^{-1/2}/ h^{3/2},\\ \frac{\alpha \gamma^{1/2}h^{3/2}}{N^{1/2}},&\;{\rm{ if }}\; \gamma^{1/2}LN^{-1/2}/h^{3/2}\le\alpha < \gamma^{1/2} LN^{3/2} /h^{3/2},\\ \frac{\alpha^2 h^3}{L N^2 },&\;{\rm{ if }}\; \alpha\ge \gamma^{1/2}L N^{3/2} /h^{3/2}. \end{cases} \end{equation} (3.75)

    Remark 3.10. It is interesting to compute the optimal delaminated length \ell in the different regimes. Collecting the results from the above computations we obtain (in the setting of (3.74) and with h\ge\gamma N^3 ) the scaling

    \begin{equation} \ell \approx \begin{cases} 0 & \;{\rm{ if }}\; \alpha < \gamma^{2/5}L^{3/5}/h,\\ \alpha^{1/3}\gamma^{-1/3}h^{4/3} N^{-1} &\;{\rm{ if }}\; \gamma^{2/5}L^{3/5}/h\le\alpha < \gamma^{1/4}h^{-1/4}N^{3/4},\\ \alpha \gamma^{-1/2} h^{3/2}N^{-3/2}&\;{\rm{ if }}\; \gamma^{1/4}h^{-1/4}N^{3/4}\le\alpha < \gamma^{1/2}h^{-3/2} L N^{3/2},\\ L ,&\;{\rm{ if }}\; \alpha\ge \gamma^{1/2}h^{-3/2} L N^{3/2} . \end{cases} \end{equation} (3.76)

    Theorem 4.1. Let \alpha\in(0, \frac\pi2] , h\le \frac14 L . Assume that u\in SBV^2_N(\omega_h; \mathbb R^2) obeys the boundary condition (2.10) and hasjump set contained in the product of an interval and a finite set,

    \begin{equation} J_u\subseteq (x_-,x_+)\times\{b_1, \dots b_{n-1}\} \end{equation} (4.1)

    (up to \mathcal H^1 -null sets) for some n\ge 1 and x_-, x_+ with -\frac L2\le x_-\le x_+\le \frac L2 .Assume additionally

    \begin{equation} \mathcal H^1(J_u) \ge c_J (n-1)(x_+-x_-) \end{equation} (4.2)

    for some c_J > 0 .Then

    \begin{equation} \mathcal {\cal E}_h[u]\ge c\min\left\{\frac{\alpha^2h^3}{L}, \frac{\alpha \gamma^{1/2} h^{3/2} }{N^{1/2}}+ \frac{\alpha^2h^3}{LN^2}, \frac{\alpha^2h^2}{N}\right\}, \end{equation} (4.3)

    where c may depend on c_J .

    As in the construction of the upper bound, n\ge1 denotes the number of layers, n-1\ge0 the number of interfaces.

    Remark 4.2. We note that the all regimes obtained here, except for the last one, are matched by upper bounds. In the proof below, we will see that this last regime corresponds to a very short delamination length. In the upper bound, this is forbidden by the injectivity constraint \ell_{\mathrm{arc}} \ge \frac{2\beta h}{n} , which however only arises from our specific construction.

    The upper bound {\cal E}_h[u] \le c\alpha^{1/3}\gamma^{2/3}h^{4/3} arises from the lower bound in (3.60). This was introduced to take into account the first term in the energy bound in (3.42). In particular, it represents the fact that choosing \ell_{\mathrm{arc}}\le\ell does not save any delamination energy. This fact is not considered in the lower bound estimate. Finally, the bound {\cal E}_h[u] \le c\alpha h^4 L^{-2}N^{-2} arises in the regime of small-angle total delamination as discussed before (3.70), where the scaling \beta\sim (\alpha h/L)^{1/2} is induced by the injectivity constraint. This is also not considered in the lower bound estimate.

    Lemma 4.3. Let \alpha\in(0, \frac\pi2] , h\le \frac14 L . Assume that u\in W^{1, 2}(\omega_h; \mathbb R^2) obeys the boundary condition (2.10). Then

    \begin{equation} {\cal E}_h[u]\ge c\frac{\alpha^2h^3}{L}. \end{equation} (4.4)

    This result follows easily from the compactness result for plates in [10]; we present here the short argument since it will be reused in the following.

    Proof. We consider the rectangles Q_k: = (hk, h(k+2))\times (0, h) , for k = 0, \dots, K: = \lfloor \frac{L}{h}-2\rfloor . By (2.9) and the Friesecke-James-Müller rigidity estimate [10], there is a universal constant c > 0 such that for every k there is R^k\in {\mathrm{SO}}(2) with

    \begin{equation} \int_{Q_k}|Du-R^k|^2 dx \le c \int_{Q_k} \operatorname{dist}^2(Du, {\mathrm{SO}}(2))dx\le c \int_{Q_k} W_{\mathrm{2D}}(Du) dx. \end{equation} (4.5)

    By (2.10), we can assume R^0 = R_{\alpha} and R^K = R_{-\alpha} . With a triangular inequality and \mathcal L^2(Q_k\cap Q_{k+1}) = h^2 we obtain

    \begin{equation} h^2 |R^k-R^{k+1}|^2\le 2 \int_{Q_k}|Du-R^k|^2 dx + 2 \int_{Q_{k+1}}|Du-R^{k+1}|^2 dx . \end{equation} (4.6)

    By a triangular inequality and Hölder's inequality,

    \begin{equation} |R^0-R^K|^2\le \left(\sum\limits_{k = 0}^{K-1} |R^k-R^{k+1}|\right)^2 \le K \sum\limits_{k = 0}^{K-1} |R^k-R^{k+1}|^2. \end{equation} (4.7)

    Combining the previous estimates and K\le L/h ,

    \begin{equation} |R^0-R^K|^2 \le c \frac{K}{h^2} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx \le c\frac{L}{h^3} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx . \end{equation} (4.8)

    With |R_\alpha-R_{-\alpha}|\ge 2\sin\alpha\ge \alpha the proof is concluded.

    Lemma 4.4. Let \alpha\in(0, \frac\pi2] , h\le \frac14 L . Assume that u\in SBV^2_N(\omega_h; \mathbb R^2) obeys the boundary condition (2.10) and hasjump set contained in the product of an interval and a finite set,

    \begin{equation} J_u\subseteq (x_-,x_+)\times\{b_1, \dots b_{n-1}\} \end{equation} (4.9)

    for some n\ge 1 and x_-, x_+ with -\frac L2\le x_- < x_+\le \frac L2 .Then, with \ell: = x_+-x_- ,

    \begin{equation} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx \ge c\min\left\{\frac{\alpha^2h^3}{L}, \frac{\alpha^2h^3}{n^2(\ell+h/n)}\right\}. \end{equation} (4.10)

    Proof. By (4.9), u\in W^{1, 2}((-L, x_-)\times(0, h)) . This set can be treated as in Lemma 4.3. We consider the rectangles Q^-_k: = (x_–(k+2)h, x_–kh)\times(0, h) , k = 0, \dots, K^-: = \lfloor \frac{x_-+L}{h}-2\rfloor , and obtain rotations R^k_-\in {\mathrm{SO}}(2) such that R_-^{K^-} = R_{\alpha} ,

    \begin{equation} |R_-^0-R_\alpha|^2\le K^-\sum\limits_{k = 0}^{K^–1} |R_-^k-R_-^{k+1}|^2\le c\frac{L}{h^3} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx. \end{equation} (4.11)

    Analogously, (4.9) also gives u\in W^{1, 2}((x_+, L)\times(0, h)) . The same computation, using the sets Q^+_k: = (x_++kh, x_++(k+2)h)\times(0, h) , leads to

    \begin{equation} |R_+^0-R_{-\alpha}|^2\le K^+\sum\limits_{k = 0}^{K^+-1} |R_+^k-R_+^{k+1}|^2\le c\frac{L}{h^3} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx . \end{equation} (4.12)

    By a triangular inequality,

    \begin{equation} |R_+^0-R_{-\alpha}|+ |R_+^0-R_-^0|+ |R_-^0-R_{\alpha}|\ge |R_\alpha-R_{-\alpha}|\ge 2\sin\alpha\ge \alpha. \end{equation} (4.13)

    If \max\{|R^0_+-R_{-\alpha}|, |R^0_–R_{\alpha}|\}\ge \frac 14\alpha , then (4.10) holds and the proof is concluded. Therefore for the rest of the proof we can assume

    \begin{equation} |R^0_–R^0_+|\ge \frac 12\alpha, \end{equation} (4.14)

    and we only need to deal with the central part of the domain. As in the proof of the upper bound we set b_0: = 0 , b_{n+1}: = h , and h_j: = b_{j+1}-b_j . We treat each layer (x_-, x_+)\times (b_j, b_{j+1}) separately, for j = 0, \dots, n . We consider the sets Q_k^j: = (x_-+(k-1)h_j, x_-+(k+1)h_j)\times(b_j, b_j+h_j) , for k = 0, \dots, K_j: = \lfloor \frac{x_+-x_-}{h_j}+1\rfloor , and obtain rotations R_j^k\in {\mathrm{SO}}(2) with

    \begin{equation} |R_j^0-R_j^{K_j}|^2\le c \frac{K_j}{h_j^2} \int_{(-L,L)\times(b_j, b_{j+1})} W_{\mathrm{2D}}(Du) dx \end{equation} (4.15)

    as well as

    \begin{equation} \int_{Q^j_0} |Du-R_j^0|^2 dx + \int_{Q^j_{K_j}} |Du-R_j^{K_j}|^2 dx \le c \int_{Q^j_0\cup Q^j_{K_j}} W_{\mathrm{2D}}(Du)dx. \end{equation} (4.16)

    Since \mathcal L^2(Q^j_0\cap Q^-_0)\ge h_j^2 and \mathcal L^2(Q^j_{K_j}\cap Q^+_0)\ge h_j^2 ,

    \begin{equation} h_j^2|R^0_j-R^0_-|^2 \le c \int_{Q^j_0} W_{\mathrm{2D}}(Du) dx + c\int_{Q^j_0\cap Q^-_0} |Du-R^0_-|^2dx, \end{equation} (4.17)

    and the same on the other side. With (4.15) we obtain

    \begin{equation} \begin{split} h_j^2|R^0_–R^0_+|^2 \le &c \int_{Q^j_0\cup Q^j_{K_j}} W_{\mathrm{2D}}(Du) dx + c\int_{Q^j_0\cap Q^-_0} |Du-R^0_-|^2dx + c\int_{Q^j_{K_j}\cap Q^+_0} |Du-R^0_+|^2dx \\ & +c K_j \int_{(-L,L)\times(b_j, b_{j+1})} W_{\mathrm{2D}}(Du) dx . \end{split} \end{equation} (4.18)

    Let A: = \{j\in\{1, \dots, n-1\}: h_j\ge \frac h{2n}\} . Then \sum_{j\not\in A} h_j\le \frac h2 , which implies

    \begin{equation} \sum\limits_{j\in A} h_j^2 \ge \frac{1}{\# A}(\sum\limits_{j\in A} h_j)^2 \ge \frac{h^2}{4n}. \end{equation} (4.19)

    We sum (4.18) over all j\in A , use that the domains of integration for different j are disjoint and that j\in A implies K_j\le 1+\ell/h_j\le 1+2n\ell/h\le 2(h+n\ell)/h , and obtain

    \begin{equation} \frac{h^2}{4n} |R^0_–R^0_+|^2 \le \sum\limits_{j\in A} h_j^2 |R^0_–R^0_+|^2 \le c (1+\frac{n\ell}{h})\int_{\omega_h} W_{\mathrm{2D}}(Du) dx. \end{equation} (4.20)

    Recalling (4.14),

    \begin{equation} \int_{\omega_h} W_{\mathrm{2D}}(Du) dx\ge c \frac{\alpha^2 h^3}{n^2(\ell+h/n)} \end{equation} (4.21)

    which concludes the proof.

    Proof of Theorem 4.1. If n = 1 or \ell = 0 then u\in W^{1, 2}(\omega_h; \mathbb R^2) and the assertion follows from Lemma 4.3. Therefore we can assume n\ge 2 , \ell > 0 .

    By (4.2) and Lemma 4.4 we have

    \begin{equation} {\cal E}_h[u]\ge c_J n\gamma \ell + c\min\left\{\frac{\alpha^2h^3}{L}, \frac{\alpha^2h^3}{n^2(\ell+h/n)}\right\}, \end{equation} (4.22)

    where \ell: = x_+-x_-\in(0, L] , n\in[1, N] . As 1/(a+b)\ge\min\{1/2a, 1/2b\} for a, b > 0 we have, with c': = \min\{c_J, c/2\} ,

    \begin{equation} {\cal E}_h[u]\ge c' \left[n\gamma \ell + \min\left\{\frac{\alpha^2h^3}{L}, \frac{\alpha^2h^3}{\ell n^2}, \frac{\alpha^2h^2}{n}\right\}\right], \end{equation} (4.23)

    which implies

    \begin{equation} {\cal E}_h[u]\ge c'\min\left\{\frac{\alpha^2h^3}{L}, n\gamma \ell + \frac{\alpha^2h^3}{\ell n^2}, \frac{\alpha^2h^2}{n}\right\}. \end{equation} (4.24)

    We treat the three terms separately. For the first one there is nothing to do, for the last one we simply use n\le N . For the middle one we use two estimates. From n\le N and \ell\le L , we get

    \begin{equation} n\gamma \ell + \frac{\alpha^2h^3}{\ell n^2}\ge \frac{\alpha^2h^3}{LN^2} \end{equation} (4.25)

    and using ab\ge2\sqrt{ab} first, and then n\le N ,

    \begin{equation} n\gamma \ell + \frac{\alpha^2h^3}{\ell n^2}\ge 2 \frac{\alpha h^{3/2}\gamma^{1/2} }{n^{1/2}}\ge 2 \frac{\alpha h^{3/2}\gamma^{1/2} }{N^{1/2}}. \end{equation} (4.26)

    Therefore

    \begin{equation} n\gamma \ell + \frac{\alpha^2h^3}{\ell n^2}\ge \max\left\{ \frac{\alpha h^{3/2}\gamma^{1/2} }{N^{1/2}}, \frac{\alpha^2h^3}{LN^2}\right\} \ge \frac12 \left( \frac{\alpha h^{3/2}\gamma^{1/2} }{N^{1/2}}+ \frac{\alpha^2h^3}{LN^2}\right), \end{equation} (4.27)

    and inserting this bound in (4.24) yields

    \begin{equation} {\cal E}_h[u]\ge \frac12 c'\min\left\{\frac{\alpha^2h^3}{L}, \frac{\alpha \gamma^{1/2} h^{3/2}}{N^{1/2}} + \frac{\alpha^2h^3}{LN^2}, \frac{\alpha^2h^2}{N}\right\} \end{equation} (4.28)

    which concludes the proof.

    Remark 4.5. We remark that the first inequality in (4.26) corresponds to the choice \ell: = \alpha h^{3/2}/(\gamma^{1/2}n^{3/2}) . Therefore in the regime in which the energy scales as \frac{\alpha \gamma^{1/2} h^{3/2} }{N^{1/2}} , delamination necessarily occurs over a length which is (up to a factor) equal to \ell: = \alpha h^{3/2}/(\gamma^{1/2}n^{3/2}) . This implies localization of delemination to a small part of the sample, in agreement with the upper bounds discussed above and the experimental observations in Figure 2.

    Analogously, the lower bound shows that bending localizes. Indeed, from (4.12) and the corresponding equation we obtain

    \begin{equation} \frac{h^3}{L}\left[ |R_+^0-R_{-\alpha}|^2+ |R_-^0-R_{\alpha}|^2\right]\le c {\cal E}_h[u]. \end{equation} (4.29)

    We recall that R_+^0-R_{-\alpha} and R_-^0-R_{\alpha} are the differences in rotation inside the two non-delaminated parts of the sample. Therefore if the energy is significantly smaller than c\alpha^2h^3/L the two non-delaminated parts of the sample carry a small part of the total rotation, and bending localizes to the delaminated part.

    We have devised a variational model for the study of the delamination of paperboard undergoing bending. As illustrated in Figure 6 (left), a rich variety of energy scaling regimes emerges (non-delaminated, locally delaminated, and totally delaminated), where the locally delaminated regime exhibits the inelastic hinge found in experiments. In the case of a small coefficient \gamma for the delamination energy, our model exhibits total delamination, where each layer deforms independently (apart from injectivity constraints). The small-angle total delamination regime is not included in the figure. This regime would replace the full delamination regime with linear growth, but only occurs for extremely small \gamma .

    The bending moment is displayed in Figure 6. We expect that the moment is continuous in the bending angle \alpha , apart from the point of the first order phase transition when delamination becomes energetically favorable. At the other scaling regime transitions, the minimizer is continuous, therefore we expect the energy to exhibit a smooth crossover with continuous moment. This figure should be compared to [4,Figure 18] – the initial undelaminated elastic response, the flat regime of increasing delamination, and the final increase of the moment are all displayed there – although the final increase seems to be of a different origin in the experiments (where it is an artefact of the set-up) than in our variational model (where it is due to total delamination, which may only occur for unphysically small \gamma ). The moment discontinuity we observe in our variational model at the onset of delamination, however, seems to be prevented by the introduction of an initial crease in the experiments in [4]. What is observed, however, is a softening behavior (see, e.g., [13], Figure 4]) in the dependence of the moment on the bending angle. A more detailed analysis of the time-dependent processes of the delamination or a cohesive zone model may regularize the discontinuity to recover a softening instead, but this is beyond the scope of this work.

    Future work should include a more detailed study of the paperboard material, taking into account the effect of its constituent cellulose fibers. This could involve also a treatment via the theory of homogenization, which was applied to the study of strength and fatigue of materials for example in [20]. From a mathematical perspective, some open problems remain for the lower bound energy estimate: currently, it is not completely ansatz-free, as we make assumptions on the delamination sets. Furthermore, the lower bound corresponding to the partially delaminated regime is not yet available.

    This work was partially supported by the Deutsche Forschungsgemeinschaft through projects 211504053/SFB1060, 441211072/SPP2256, and 441523275/SPP2256. PD and JO are grateful for the hospitality afforded by the Institute for Applied Mathematics of the University of Bonn.

    The authors declare no conflict of interest.



    [1] L. Ambrosio, N. Fusco, D. Pallara, Functions of bounded variation and free discontinuity problems, New York: The Clarendon Press, Oxford University Press, 2000.
    [2] S. S. Antman, Nonlinear problems of elasticity, New York: Springer, 2005. https://doi.org/10.1007/0-387-27649-1
    [3] B. Audoly, Y. Pomeau, Elasticity and geometry: From hair curls to the non-linear response of shells, Oxford: Oxford University Press, 2010.
    [4] L. A. A. Beex, R. H. J. Peerlings, An experimental and computational study of laminated paperboard creasing and folding, Int. J. Solids Struct., 46 (2009), 4192–4207. https://doi.org/10.1016/j.ijsolstr.2009.08.012 doi: 10.1016/j.ijsolstr.2009.08.012
    [5] L. A. A. Beex, R. H. J. Peerlings, On the influence of delamination on laminated paperboard creasing and folding, Phil. Trans. R. Soc. A, 370 (2012), 1912–1924. https://doi.org/10.1098/rsta.2011.0408 doi: 10.1098/rsta.2011.0408
    [6] H. Ben Belgacem, S. Conti, A. DeSimone, S. Müller, Rigorous bounds for the Föppl-von Kármán theory of isotropically compressed plates, J. Nonlinear Sci., 10 (2000), 661–683. https://doi.org/10.1007/s003320010007 doi: 10.1007/s003320010007
    [7] B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, J. Elasticity, 91 (2008), 5–148. https://doi.org/10.1007/s10659-007-9107-3 doi: 10.1007/s10659-007-9107-3
    [8] D. P. Bourne, S. Conti, S. Müller, Energy bounds for a compressed elastic film on a substrate, J. Nonlinear Sci., 27 (2017), 453–494. https://doi.org/10.1007/s00332-016-9339-0 doi: 10.1007/s00332-016-9339-0
    [9] S. Conti, F. Maggi, Confining thin elastic sheets and folding paper, Arch. Rational Mech. Anal., 187 (2008), 1–48. https://doi.org/10.1007/s00205-007-0076-2 doi: 10.1007/s00205-007-0076-2
    [10] G. Friesecke, R. D. James, S. Müller, A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity, Commun. Pure Appl. Math., 55 (2002), 1461–1506. https://doi.org/10.1002/cpa.10048 doi: 10.1002/cpa.10048
    [11] H. Huang, Numerical and experimental investigation of paperboard creasing and folding, PhD thesis, KTH Royal Institute of Technology, 2011.
    [12] W. Jin, P. Sternberg, Energy estimates of the von Kármán model of thin-film blistering, J. Math. Phys., 42 (2001), 192–199. https://doi.org/10.1063/1.1316058 doi: 10.1063/1.1316058
    [13] M. Klingenberg, A. Boldizar, K. Hofer, Mechanical properties of paperboard with a needled middle layer, Cell. Chem. Technol., 52 (2018), 89–97.
    [14] R. V. Kohn, S. Müller, Surface energy and microstructure in coherent phase transitions, Commun. Pure Appl. Math., 47 (1994), 405–435. https://doi.org/10.1002/cpa.3160470402 doi: 10.1002/cpa.3160470402
    [15] B. Lawn, Fracture of brittle solids, Cambridge: Cambridge University Press, 1993. https://doi.org/10.1017/CBO9780511623127
    [16] M. Lecumberry, S. Müller, Stability of slender bodies under compression and validity of the von Kármán theory, Arch. Rational Mech. Anal., 193 (2009), 255–310. https://doi.org/10.1007/s00205-009-0232-y doi: 10.1007/s00205-009-0232-y
    [17] H. LeDret, A. Raoult, The nonlinear membrane model as a variational limit of nonlinear three-dimensional elasticity, J. Math. Pure. Appl., 73 (1995), 549–578.
    [18] E. Linvill, M. Wallmeier, S. Östlund, A constitutive model for paperboard including wrinkle prediction and post-wrinkle behavior applied to deep drawing, Int. J. Solids Struct., 117 (2017), 143–158. https://doi.org/10.1016/j.ijsolstr.2017.03.029 doi: 10.1016/j.ijsolstr.2017.03.029
    [19] A. E. Lobkovsky, S. Gentges, H. Li, D. Morse, T. A. Witten, Scaling properties of stretching ridges in a crumpled elastic sheet, Science, 270 (1995), 1482–1485. https://doi.org/10.1126/science.270.5241.1482 doi: 10.1126/science.270.5241.1482
    [20] J. Orlik, Homogenization of strength, fatigue and creep durability of composites with near periodic structure, Math. Mod. Meth. Appl. Sci., 15 (2005), 1329–1347. https://doi.org/10.1142/S0218202505000807 doi: 10.1142/S0218202505000807
    [21] S. Östlund, Three-dimensional deformation and damage mechanisms in forming of advanced structures in paper, In: Advances in pulp and paper research, Trans. of the XVIth Fund. Res. Symp. Oxford, 2017,489–594.
    [22] J.-W. Simon, A review of recent trends and challenges in computational modeling of paper and paperboard at different scales, Arch. Computat. Methods Eng., 28 (2021), 2409–2428. https://doi.org/10.1007/s11831-020-09460-y doi: 10.1007/s11831-020-09460-y
    [23] S. C. Venkataramani, Lower bounds for the energy in a crumpled elastic sheet – a minimal ridge, Nonlinearity, 17 (2004), 301–312. https://doi.org/10.1088/0951-7715/17/1/017 doi: 10.1088/0951-7715/17/1/017
  • This article has been cited by:

    1. R. Cavuoto, A. Cutolo, K. Dayal, L. Deseri, M. Fraldi, Distal and non-symmetrical crack nucleation in delamination of plates via dimensionally-reduced peridynamics, 2023, 172, 00225096, 105189, 10.1016/j.jmps.2022.105189
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1952) PDF downloads(129) Cited by(1)

Figures and Tables

Figures(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog