
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
[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] | 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 |
[5] | 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 |
[6] | 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 |
[7] | 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 |
[8] | 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 |
[9] | 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 |
[10] | Helen Christodoulidi, Christos Efthymiopoulos . Stages of dynamics in the Fermi-Pasta-Ulam system as probed by the first Toda integral. Mathematics in Engineering, 2019, 1(2): 359-377. doi: 10.3934/mine.2019.2.359 |
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.
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 N−1 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:Ωh→R3 which jump only on (a subset of) N−1 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:x2∈hNZ}=(−L,L)×{hN,2hN,…,(N−1)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 ∇U∈L2(Ω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:Ωh→R3 such that the distributional gradient DU is a measure of the form DU=∇UL3+[U]⊗νH2∟JU, with ∇U∈L1(Ωh;R3×3), JU the H2-rectifiable jump set of U, [U]:JU→R3 its jump, and ν:JU→S2 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 U∈SBV2N(Ω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),forU∈SBV2N(Ω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 u∈SBV2N(ωh;R2), with ωh:=(−L,L)×(0,h). The latter set is defined as the set of SBV functions such that ∇u∈L2(ωh;R2×2) and
Ju⊂{x∈ωh:x2∈hNZ}=(−L,L)×{hN,2hN,…,(N−1)hN}, | (2.7) |
and one easily sees that U∈SBV2N(Ωh;R3) if and only if u∈SBV2N(ωh;R2).
The energy then reduces to
Eh[u]:=E3Dh[U]=∫(−L,L)×(0,h)W2D(∇u)dx+γH1(Ju),foru∈SBV2N(ωh;R2), | (2.8) |
where W2D is defined by W2D(ξ):=W3D(ξ+e3⊗e3) 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 x2−h2). 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 u∈W1,∞(ω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),if−L<x1<−L2,f(x1),if−L2≤x1≤L2,f(L2)+(x1−L2)f′(L2),ifL2<x1<L. | (3.3) |
One easily verifies that v∈W2,∞((−L,L);R2), with |v′|=1 and |v"|≤2α/L almost everywhere. We then define u:ωh→R2 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 u∈W1,∞(ω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]≤∫L−L∫h0c|x2v"(x1)|2dx2dx1=L∫h0c(2αx2L)2dx2≤cα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 {x1≥0}, provided that we impose the condition ucpa1(0,x2)=0. We fix a parameter ζ∈(0,L/(2h)], and assume that
Ducpa=R−αforζ(h−x2)<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<ζ(h−x2) 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 a1∈R (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 F∗12=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=ζ(h−x2)}, 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α)=1−cosα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α(β):=1−cosα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, βeq∈C1((0,π2)), and from
∂fα(β)∂α=sinβcosβ1−cos(α−β)>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
1−cosα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:¯ωh→R2 defined by
ucpa(x):={(x1cosβ|x1|sinβ+dx2),if|x1|<ζ(h−x2),(x1cosα+(x2−h)sinαsgnx1−|x1|sinα+(x2−h)cosα+dh),ifζ(h−x2)≤|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
(1−cosα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)=(1−cosβ)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<ζ(h−x2), Ducpa=R−α for x1≥ζ(h−x2), 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.
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 n≤N layers of paperboard which have been separated across n−1 delamination surfaces. We assume n≥1; the case n=1 without delamination has already been treated in Section 3.1. To this end, we fix a subset {b1<b2<⋯<bn−1}⊆(0,h)∩hNZ and assume that delamination occurs only on the surfaces {x2=b1,…bn−1}, in the sense that Ju⊆(−L,L)×{b1,…,bn−1}. For notational convenience we denote b0:=0 and bn:=h. We label by hj:=bj+1−bj, 0≤j<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,…,ˆfn−1:[−L,L]→R2 by
ˆfj(x1):=ucpa(x1,bj)={(x1cosβ|x1|sinβ+dbj),if|x1|<lj,(x1cosα+(bj−h)sinαsgnx1−|x1|sinα+(bj−h)cosα+dh),iflj≤|x1|≤L, | (3.25) |
where lj:=ζ(h−bj) 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 ˆf′j, 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 f′j(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 ℓarc→0, 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 2ℓarc+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 φj∈W1,∞(R) by
φj(x1):={x1ℓarcβ,if|x1|≤ℓarc,βsgnx1,ifℓarc<|x1|<lj+ℓarc,βsgnx1−|x1|−lj−ℓarcℓarc(α+β)sgnx1,iflj+ℓarc≤|x1|≤lj+2ℓarc,−αsgnx1,if|x1|>lj+2ℓarc, | (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(x′1)e1dx′1. | (3.27) |
We observe that fj∈W2,∞([−L,L];R2), with |f′j|=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), f′j(x1+ℓarc)=ˆf′j(x1)=Rβe1 for x1∈(0,lj) and f′j(x1+2ℓarc)=ˆf′j(x1)=R−αe1 for x1>lj. A short computation shows that if x1≥max{lj+2ℓarc,lj+1+2ℓarc} then
fj+1(x1)−fj(x1)=ˆfj+1(0)−ˆfj(0)+(lj+1−lj)(Rβe1−R−αe1)=ˆfj+1(x1)−ˆfj(x1)=R−αe2hj. | (3.28) |
Further, we observe that
lj≤l0=hζ=hsinαcosα−cosβ. | (3.29) |
Step 3. We are now ready to define the required mapping u:ωh→R2. 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)+(x2−bj)(f′j)⊥(x1),forx2∈[bj,bj+hj). | (3.30) |
This map clearly belongs to SBV2N(ωh;R2).
Assume now that l0+2ℓarc≤12L. Then for x1≥L2 we have f′j=R−αe1, and (3.28) holds. Therefore u is continuous for x1≥L2, 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]≤n−1∑j=12γ(lj+2ℓarc)+n−1∑j=0∫L−L∫bj+1bjc|(x2−bj)fj"(x1)|2dx2dx1≤4nγ(hζ+ℓarc)+cn−1∑j=0β2h3jℓarc. | (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
hj≤ℓarcβforj=0,…,n−1. | (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)+λ(f′j)⊥(x1)≠fj+1(x′1)foralljandx1,x′1∈R,λ∈[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)2≥cos2βsin2βA22+(A22+2A2dhj+d2h2j)=1sin2βA22+2A2dhj+d2h2j≥d2h2j−d2h2jsin2β=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.
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 δ:=tan−1(ζ) and orthogonal unit vectors ˉe1:=Rα−δe1, ˉe2:=Rα−δe2=(sin(δ−α)cos(δ−α)). From (3.28) and lj=lj+1+ζhj we have, for t+ζhj≤−L/2, fj+1(t+ζhj)−fj(t)=hjRα(e2+ζe1)=hj√1+ζ2ˉe2. As f′j+1(t+ζhj)=f′j(t) for t+ζhj<−ℓarc, we conclude that in this range fj+1(t+ζhj)=fj(t)+hj√1+ζ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δ=hj√1+ζ2hj. The angle ˉδ between the down-slope and ˉe1 satisfies cosˉδ=(dcosβ)hj√1+ζ2hj≥cosδ, since β≤βeq and therefore dcosβ≥1. For an illustration, see Figure 5a. Alternatively, (3.37) can be proven algebraically. First, φj∈[−β,α] for x1≤0, 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(1tanm−1tanp). | (3.38) |
On the other hand, by (3.16) the assumption β≤βeq (in the form fα(β)≥1) is the same as
1tanm≥1+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,t≤0. Setting ˉA:=(ˉe1⋅(fj(s)−fj(t))ˉe2⋅(fj(s)−fj(t))), we calculate
|fj+1(t+hjζ)−fj(s)|2=ˉA21+(ˉA2+√1+ζ2hj)2≥(1tan2δ+1)ˉA22+2ˉA2√1+ζ2hj+(1+ζ2)h2j=1+ζ2ζ2ˉA22+2ˉA2√1+ζ2hj+(1+ζ2)h2j≥h2j. | (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,√1+ζ2h}.
Lemma 3.3. Fix α∈(0,π4], h>0, L≥h, N∈N, N≥1. For any n∈N with 1≤n≤N, any β∈(α,βeq(α)], and any ℓarc>0, if
2βhn≤ℓarc≤18Landℓ:=sinαcosα−cosβh≤14L | (3.41) |
then there exists a map u∈SBV2N(ωh;R2), which obeys (2.10), is injective, and such that
Eh[u]≤c(γ(ℓ+ℓarc)n+β2h3ℓarcn2)=c(γn(hsinαcosα−cosβ+ℓarc)+β2h3ℓarcn2) | (3.42) |
with a constant c>0 only depending on the elastic energy density W2D.
Proof. For n=1 there is no delamination, the assertion follows from Lemma 3.1, α≤β, and ℓarc≤L. Therefore we can assume n≥2 in the following.
We first choose a subset {b1,…,bn−1}⊆{h1N,…,hN−1N} such that, setting b0=0 and bn=h, we have hj=bj+1−bj≤2hn for all j=0,…,n−1. We immediately notice that the map u given in (3.30) is a member of SBV2N(ωh;R2), with its jump set contained in the set (−2ℓarc−ℓ,2ℓarc+ℓ)×{b1,…,bn−1}. The energy estimate follows then from (3.31).
The conditions in (3.41) imply in particular that 2ℓarc+ℓ≤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 hj≤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 η∈(0,1] such that for 0<h≤ηL, α∈(0,π/2] there is u∈SBV2N(Ω;R2) which obeys (2.10), u∈W1,∞(Ω∖¯Ju;R2) with H1(¯Ju∖Ju)=0, is injective, and obeys
Eh[u]≤Cmin{α2h3L,α1/3γ2/3h4/3+αγ1/2h3/2N1/2+α2h3LN2+αh4L2N2}. | (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 α is large, the second treats the case of small α. As h3LN2≤h4L2N2, in the case of large α the last term does not appear.
Proposition 3.6. There are C>0 and η∈(0,1) such that for 0<h≤ηL, α∈(0,π/2] there is u∈SBV2N(Ω;R2) which obeys (2.10), is W1,∞ in Ω∖¯Ju, injective, and obeys
Eh[u]≤Cmin{h3L,γ2/3h4/3+γ1/2h3/2N1/2+h3LN2}. | (3.44) |
Proof. We prove the proposition in 3 steps.
Step 1. If α>π4 we can combine two folds with L′:=L/4 and α′:=α/2. Therefore we assume in the remainder of the proof that α∈(0,π4].
Step 2. By Lemma 3.1 there is an injective map u∈W1,∞(ωh;R2) with Eh[u]≤Ch3/L. This proves the first bound.
Step 3. We intend to use Lemma 3.3 with β:=βeq(α). We first check the second condition in (3.41). The function
α↦sinαcosα−cosβeq(α) | (3.45) |
is continuous on (0,π4], and as βeq=O(α1/3) for small α, it converges to 0 for α→0. Therefore there is C>0 such that
0<sinαcosα−cosβeq(α)≤Cforallα∈(0,π4]. | (3.46) |
By Lemma 3.3 we thus know that for any n∈{1,…,N} and ℓarc>0, if
2βhn≤ℓarc≤18LandCh≤14L | (3.47) |
then
infuEh[u]≤c(γhn+γℓarcn+h3ℓarcn2). | (3.48) |
We assume that
hL≤min{14C,18π}, | (3.49) |
which is guaranteed if η is chosen appropriately. Noting that reducing ℓarc below O(h) does not reduce the energy, due to the first term in (3.48), we can choose
ℓarc∈[πh,L8]. | (3.50) |
Since β≤π2, we are automatically assured that all requirements in (3.47) are satisfied. At the same time, if ℓ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 ℓarc is nonempty.
We first optimize the choice of n. Set
n:=min{N,⌈hγ1/3ℓ2/3arc⌉}∈{1,…,N}. | (3.51) |
Inserting in (3.48) leads to
infuEh[u]≤{c(γℓarc+γ2/3hℓ1/3arc),ifh≤Nγ1/3ℓ2/3arc,ch3ℓarcN2,otherwise, | (3.52) |
where the term γℓarc stems from the fact the n≥1.
It remains to choose the value of ℓ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 ℓarccrit:=γ−1/2h3/2N−3/2∈[πh,L/8] then we obtain a bound of the form h3/2γ1/2/N1/2. Otherwise, the relevant constraint is the one that prevents this optimal value from being chosen. Specifically,
infuEh[u]≤{c(γh+γ2/3h4/3),ifγ−1/2h3/2N−3/2<πh,cγ1/2h3/2N1/2,ifγ−1/2h3/2N−3/2∈[πh,L/8],ch3LN2,ifL<8γ−1/2h3/2N−3/2, | (3.53) |
or, equivalently,
infuEh[u]≤{c(γh+γ2/3h4/3),ifh<γN3,cγ1/2h3/2N1/2,ifγN3≤h≤γ1/3L2/3N,ch3LN2,ifγ1/3L2/3N<h. | (3.54) |
We remark that the first and the last regime are disjoint, since γ1/3L2/3N<h and h≤L imply γ1/3h2/3N<h, which is equivalent to γN3<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 γh+γ2/3h4/3 has the shortest possible delamination length ℓarc=πh and n∼⌈γ−1/3h1/3⌉. It originates by balancing the two terms γℓarcn+h3ℓarcn2 after setting ℓarc=h.
The term γh, corresponding to n=1, can be dropped. To see this, we first note that it is only relevant if h≤γ. In this regime, however, it would be convenient not to delaminate at all, obtaining an energy h3L−1. Indeed, as h≤L, h≤γ implies h3L−1≤h2≤γh.
Localized full delamination: The energy scales as γ1/2h3/2N1/2. As n=N, each layer is delaminated, however only over a length ℓarc=γ−1/2h3/2N−3/2. This originates by balancing the two terms γℓarcn+h3ℓarcn2 after setting n=N. All layers are delaminated, but only over a length ℓarc. The energy corresponds to N plates of thickness h/N bent over a length ℓarc, the value is determined so that this exactly balances the delamination energy.
Total delamination: The energy scales as h3LN2, 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, α0∈(0,π/2] such that for 0<h≤L, α∈(0,α0] there is u∈SBV2N(Ω;R2) which obeys (2.10), is W1,∞ in Ω∖¯Ju, injective, and obeys
Eh[u]≤Cmin{α2h3L,α1/3γ2/3h4/3+αγ1/2h3/2N1/2+α2h3LN2+αh4L2N2}. | (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∈W1,∞(ωh;R2) with Eh[u]≤Cα2h3/L.
Step 2. The delaminated construction is obtained with Lemma 3.3 with a careful choice of the parameters. In this regime of small α, any delaminated construction must also satisfy that β is small, as β≤βeq(α)≃(4α)1/3 for α→0. A straightforward Taylor series expansion of the delamination length factor ζ leads to
ζ=sinαcosα−cosβ=2α+O(α3)β2−α2+O(β4)+O(α4). | (3.56) |
Since α≤β, the O(α4) in the denominator can be ignored. It is however important to avoid that β2−α2 in the denominator becomes too small. For this reason we shall restrict the choice of β by assuming β≥2α, which implies cosα−cosβ≥cosβ2−cosβ=38β2+O(β4). For β sufficiently small, this is larger than β2/3. Therefore there is α0>0 such that for all α∈(0,α0] one has 2α≤βeq(α) and for all β∈[2α,βeq(α)]
ζ=sinαcosα−cosβ≤3αβ2. | (3.57) |
By Lemma 3.3 we know that for any n∈{1,…,N} and ℓarc>0, if
2βhn≤ℓarc≤18Land3αβ2h≤14L | (3.58) |
(noting that h≤L guarantees that the set of admissible β≤βeq is not empty, after potentially decreasing α0 again) then
infuEh[u]≤c(αβ2γhn+γℓarcn+β2h3ℓarcn2). | (3.59) |
We assume that
αβ2h≤ℓarc≤18L | (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 β≤βeq(α)≤cα1/3, and n≥1, we obtain for some c∗>0
2βhn≤c∗αβ2h. | (3.61) |
Without loss of generality we can assume c∗≥32. We then see that all conditions in (3.58) and in (3.60) are satisfied provided
c∗αβ2h≤ℓarc≤18L. | (3.62) |
By the same argument we see that this condition is optimal (up to a factor). We choose
n:=min{N,⌈β2/3hγ1/3ℓ2/3arc⌉}∈{1,…,N}. | (3.63) |
Inserting in (3.59) leads to
infuEh[u]≤{c(γℓarc+β2/3γ2/3hℓ1/3arc),ifβ2/3h≤γ1/3ℓ2/3arcN,cβ2h3ℓarcN2,otherwise, | (3.64) |
noting again that the term γℓarc stems from rounding n up to the next integer, which is needed as we require n≥1. As in Step 3 of the proof of Proposition 3.6, the first expression is increasing in ℓarc, the second decreasing. Therefore the optimal value is the critical one, ℓarccrit:=βγ−1/2h3/2N−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 ℓarc=ℓarccrit is admissible. Then one obtains n=N, and infuEh[u]≤cβγ1/2h3/2N1/2. Restating the requirements on ℓarc from (3.62), this is possible if ℓarccrit≥c∗αh/β2 and ℓarccrit≤L/8, which is the same as
h≥c2∗α2β−6γN3 | (3.65) |
and
h≤14β−2/3γ1/3L2/3N. | (3.66) |
Assume now that one of these two conditions is violated (one can easily check that if c∗αβ2h≤18L, so that (3.62) gives a nonempty set of admissible values of ℓ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 ℓarc=c∗αβ−2h, n=⌈β2/3hγ1/3ℓ2/3arc⌉, and infuEh[u]≤c(αβ−2γh+α1/3γ2/3h4/3).
● If (3.66) is violated, then ℓarc=L/8, n=N, and infuEh[u]≤cβ2h3LN2.
Summarizing, this leads to
E:=infuEh[u]≤{c(αγhβ2+α1/3γ2/3h4/3),ifh<c2∗α2β−6γN3,cβγ1/2h3/2N1/2,ifc2∗α2β−6γN3≤h≤14β−2/3γ1/3L2/3N,cβ2h3LN2,ifh>14β−2/3γ1/3L2/3N, | (3.67) |
or, equivalently,
E≤{c(αγhβ2+α1/3γ2/3h4/3),ifβ<α1/3γ1/6h−1/6N1/2,cβγ1/2h3/2N1/2,ifα1/3γ1/6h−1/6N1/2≤β≤γ1/2h−3/2LN3/2,cβ2h3LN2,ifβ>γ1/2h−3/2LN3/2. | (3.68) |
It remains to choose β∈[2α,βeq(α)] with c∗αβ−2h≤18L (from (3.62)), i.e.,
max{2α,(8c∗αh/L)1/2}≤β≤βeq(α). | (3.69) |
These bounds scale as α+α1/2(h/L)1/2 and α1/3, respectively, and in particular, possibly reducing α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 βcrit:=α1/3γ1/6h−1/6N1/2, and that – provided α0 is chosen sufficiently small – α1/3≤βeq(α) for all α∈(0,α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 β=βeq(α), 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 βcrit>α1/3, which is the same as h<γN3, then we take β=α1/3, and we are in the first regime in (3.68). This gives E≤cα1/3γh+cα1/3γ2/3h4/3. The first term does not contribute. Indeed, if α1/3γ2/3h4/3≤α2h3/L then necessarily γ2/3L≤α5/3h5/3, which, since h≤L and α≤1, implies γ≤h. Therefore we can drop the α1/3γh term and obtain E≤cmin{α1/3γ2/3h4/3,α2h3/L}. In this regime ℓarc=α1/3h and n≃γ−1/3h1/3≤N, so that the energy bound behaves as γℓarcn.
Sharp fold full delamination: If max{2α,(8c∗αh/L)1/2}≤βcrit≤α1/3, then β=βcrit is optimal and E≤cα1/3γ2/3h4/3. In this regime n=N and ℓarc=α1/3γ−1/3h4/3N−1 (indeed, the bound on E is of the form γℓarcN).
If βcrit<max{2α,(8c∗αh/L)1/2} then we take β=max{2α,(8c∗α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 α≥2c∗h/L, so that β=2α. If 2α≤γ1/2h−3/2LN3/2 then we are in the second regime in (3.68) (with n=N and ℓarc=ℓarccrit=2αγ−1/2h3/2N−3/2) and E≤cαγ1/2h3/2N1/2.
Total delamination: If α≥2c∗h/L and 2α>γ1/2h−3/2LN3/2, then we are in the third regime in (3.68) (with n=N and ℓarc=L/8) and E≤cα2h3LN2.
Small-angle total delamination: If instead α<2c∗h/L then β=(8c∗αh/L)1/2. A short computation shows that βcrit≤(αh/L)1/2 is equivalent to (αh/L)1/2≥γ1/2h−3/2LN3/2, so that we are necessarily in the third regime in (3.68) (with n=N and ℓarc=L/8) and E≤cαh4L2N2.
In all cases we may conclude
infuEh[u]≤cmin{α2h3L,α1/3γ2/3h4/3+αγ1/2h3/2N1/2+αh4L2N2+α2h3LN2}. | (3.70) |
Remark 3.8. A more precise summary of the cases above shows that for h≤γN3, then only the "localized partial delamination" regime is relevant, and one obtains (for α∈(0,α0])
infuEh[u]≤cmin{α2h3L,α1/3γ2/3h4/3}. | (3.71) |
Recalling (3.54) one can see that this holds for all α∈(0,π2]. As the delamination length is proportional to α1/3h, it is always smaller than L. In turn, the number of delaminated layers n≃γ−1/3h1/3 does not depend on the opening α. Both the delamination length and the number of delaminated layers are discontinuous at yielding point α≃γ2/5L3/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
Eh[u]≤Cmin{α2h3L,max{α1/3γ2/3h4/3,αγ1/2h3/2N1/2,αh4L2N2,α2h3LN2}}. | (3.72) |
To see this, it suffices to observe that max{a,b}∼a+b for any a,b≥0.
In order to gain some understanding in these many regimes we sort them in increasing values of α (see also Figure 6 below). For definiteness let us assume that the first linear term is the relevant one, in the sense that γ1/2h3/2N−1/2≥h4L−2N−2, so that the last term in the maximum can be ignored. The condition is equivalent to h5≤γL4N3. A simple computation shows that the maximum is then equal to
{α1/3γ2/3h4/3,ifα<γ1/4N3/4/h1/4,αγ1/2h3/2N1/2,ifγ1/4N3/4/h1/4≤α≤γ1/2LN3/2/h3/2,α2h3LN2,ifα>γ1/2LN3/2/h3/2. | (3.73) |
It remains to take the minimum between this expression and α2h3/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 γ2/5L3/5/h<γ1/4N3/4/h1/4 (which is the same as γL4<h5N5) then we obtain
1cEh[u]≤{α2h3L,ifα<γ2/5L3/5/h,α1/3γ2/3h4/3,ifγ2/5L3/5/h≤α<γ1/4N3/4/h1/4,αγ1/2h3/2N1/2,ifγ1/4N3/4/h1/4≤α<γ1/2LN3/2/h3/2,α2h3LN2,ifα≥γ1/2LN3/2/h3/2. | (3.74) |
whereas for γ2/5L3/5/h≥γ1/4N3/4/h1/4 the α1/3 regime disappears and
1cEh[u]≤{α2h3L,ifα<γ1/2LN−1/2/h3/2,αγ1/2h3/2N1/2,ifγ1/2LN−1/2/h3/2≤α<γ1/2LN3/2/h3/2,α2h3LN2,ifα≥γ1/2LN3/2/h3/2. | (3.75) |
Remark 3.10. It is interesting to compute the optimal delaminated length ℓ in the different regimes. Collecting the results from the above computations we obtain (in the setting of (3.74) and with h≥γN3) the scaling
ℓ≈{0ifα<γ2/5L3/5/h,α1/3γ−1/3h4/3N−1ifγ2/5L3/5/h≤α<γ1/4h−1/4N3/4,αγ−1/2h3/2N−3/2ifγ1/4h−1/4N3/4≤α<γ1/2h−3/2LN3/2,L,ifα≥γ1/2h−3/2LN3/2. | (3.76) |
Theorem 4.1. Let α∈(0,π2], h≤14L. Assume that u∈SBV2N(ωh;R2) obeys the boundary condition (2.10) and hasjump set contained in the product of an interval and a finite set,
Ju⊆(x−,x+)×{b1,…bn−1} | (4.1) |
(up to H1-null sets) for some n≥1 and x−,x+ with −L2≤x−≤x+≤L2.Assume additionally
H1(Ju)≥cJ(n−1)(x+−x−) | (4.2) |
for some cJ>0.Then
Eh[u]≥cmin{α2h3L,αγ1/2h3/2N1/2+α2h3LN2,α2h2N}, | (4.3) |
where c may depend on cJ.
As in the construction of the upper bound, n≥1 denotes the number of layers, n−1≥0 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 ℓarc≥2βhn, which however only arises from our specific construction.
The upper bound Eh[u]≤cα1/3γ2/3h4/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 ℓarc≤ℓ does not save any delamination energy. This fact is not considered in the lower bound estimate. Finally, the bound Eh[u]≤cαh4L−2N−2 arises in the regime of small-angle total delamination as discussed before (3.70), where the scaling β∼(αh/L)1/2 is induced by the injectivity constraint. This is also not considered in the lower bound estimate.
Lemma 4.3. Let α∈(0,π2], h≤14L. Assume that u∈W1,2(ωh;R2) obeys the boundary condition (2.10). Then
Eh[u]≥cα2h3L. | (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 Qk:=(hk,h(k+2))×(0,h), for k=0,…,K:=⌊Lh−2⌋. 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 Rk∈SO(2) with
∫Qk|Du−Rk|2dx≤c∫Qkdist2(Du,SO(2))dx≤c∫QkW2D(Du)dx. | (4.5) |
By (2.10), we can assume R0=Rα and RK=R−α. With a triangular inequality and L2(Qk∩Qk+1)=h2 we obtain
h2|Rk−Rk+1|2≤2∫Qk|Du−Rk|2dx+2∫Qk+1|Du−Rk+1|2dx. | (4.6) |
By a triangular inequality and Hölder's inequality,
|R0−RK|2≤(K−1∑k=0|Rk−Rk+1|)2≤KK−1∑k=0|Rk−Rk+1|2. | (4.7) |
Combining the previous estimates and K≤L/h,
|R0−RK|2≤cKh2∫ωhW2D(Du)dx≤cLh3∫ωhW2D(Du)dx. | (4.8) |
With |Rα−R−α|≥2sinα≥α the proof is concluded.
Lemma 4.4. Let α∈(0,π2], h≤14L. Assume that u∈SBV2N(ωh;R2) obeys the boundary condition (2.10) and hasjump set contained in the product of an interval and a finite set,
Ju⊆(x−,x+)×{b1,…bn−1} | (4.9) |
for some n≥1 and x−,x+ with −L2≤x−<x+≤L2.Then, with ℓ:=x+−x−,
∫ωhW2D(Du)dx≥cmin{α2h3L,α2h3n2(ℓ+h/n)}. | (4.10) |
Proof. By (4.9), u∈W1,2((−L,x−)×(0,h)). This set can be treated as in Lemma 4.3. We consider the rectangles Q−k:=(x–(k+2)h,x–kh)×(0,h), k=0,…,K−:=⌊x−+Lh−2⌋, and obtain rotations Rk−∈SO(2) such that RK−−=Rα,
|R0−−Rα|2≤K−K–1∑k=0|Rk−−Rk+1−|2≤cLh3∫ωhW2D(Du)dx. | (4.11) |
Analogously, (4.9) also gives u∈W1,2((x+,L)×(0,h)). The same computation, using the sets Q+k:=(x++kh,x++(k+2)h)×(0,h), leads to
|R0+−R−α|2≤K+K+−1∑k=0|Rk+−Rk+1+|2≤cLh3∫ωhW2D(Du)dx. | (4.12) |
By a triangular inequality,
|R0+−R−α|+|R0+−R0−|+|R0−−Rα|≥|Rα−R−α|≥2sinα≥α. | (4.13) |
If max{|R0+−R−α|,|R0–Rα|}≥14α, then (4.10) holds and the proof is concluded. Therefore for the rest of the proof we can assume
|R0–R0+|≥12α, | (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 b0:=0, bn+1:=h, and hj:=bj+1−bj. We treat each layer (x−,x+)×(bj,bj+1) separately, for j=0,…,n. We consider the sets Qjk:=(x−+(k−1)hj,x−+(k+1)hj)×(bj,bj+hj), for k=0,…,Kj:=⌊x+−x−hj+1⌋, and obtain rotations Rkj∈SO(2) with
|R0j−RKjj|2≤cKjh2j∫(−L,L)×(bj,bj+1)W2D(Du)dx | (4.15) |
as well as
∫Qj0|Du−R0j|2dx+∫QjKj|Du−RKjj|2dx≤c∫Qj0∪QjKjW2D(Du)dx. | (4.16) |
Since L2(Qj0∩Q−0)≥h2j and L2(QjKj∩Q+0)≥h2j,
h2j|R0j−R0−|2≤c∫Qj0W2D(Du)dx+c∫Qj0∩Q−0|Du−R0−|2dx, | (4.17) |
and the same on the other side. With (4.15) we obtain
h2j|R0–R0+|2≤c∫Qj0∪QjKjW2D(Du)dx+c∫Qj0∩Q−0|Du−R0−|2dx+c∫QjKj∩Q+0|Du−R0+|2dx+cKj∫(−L,L)×(bj,bj+1)W2D(Du)dx. | (4.18) |
Let A:={j∈{1,…,n−1}:hj≥h2n}. Then ∑j∉Ahj≤h2, which implies
∑j∈Ah2j≥1#A(∑j∈Ahj)2≥h24n. | (4.19) |
We sum (4.18) over all j∈A, use that the domains of integration for different j are disjoint and that j∈A implies Kj≤1+ℓ/hj≤1+2nℓ/h≤2(h+nℓ)/h, and obtain
h24n|R0–R0+|2≤∑j∈Ah2j|R0–R0+|2≤c(1+nℓh)∫ωhW2D(Du)dx. | (4.20) |
Recalling (4.14),
∫ωhW2D(Du)dx≥cα2h3n2(ℓ+h/n) | (4.21) |
which concludes the proof.
Proof of Theorem 4.1. If n=1 or ℓ=0 then u∈W1,2(ωh;R2) and the assertion follows from Lemma 4.3. Therefore we can assume n≥2, ℓ>0.
By (4.2) and Lemma 4.4 we have
Eh[u]≥cJnγℓ+cmin{α2h3L,α2h3n2(ℓ+h/n)}, | (4.22) |
where ℓ:=x+−x−∈(0,L], n∈[1,N]. As 1/(a+b)≥min{1/2a,1/2b} for a,b>0 we have, with c′:=min{cJ,c/2},
Eh[u]≥c′[nγℓ+min{α2h3L,α2h3ℓn2,α2h2n}], | (4.23) |
which implies
Eh[u]≥c′min{α2h3L,nγℓ+α2h3ℓn2,α2h2n}. | (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≤N. For the middle one we use two estimates. From n≤N and ℓ≤L, we get
nγℓ+α2h3ℓn2≥α2h3LN2 | (4.25) |
and using ab≥2√ab first, and then n≤N,
nγℓ+α2h3ℓn2≥2αh3/2γ1/2n1/2≥2αh3/2γ1/2N1/2. | (4.26) |
Therefore
nγℓ+α2h3ℓn2≥max{αh3/2γ1/2N1/2,α2h3LN2}≥12(αh3/2γ1/2N1/2+α2h3LN2), | (4.27) |
and inserting this bound in (4.24) yields
Eh[u]≥12c′min{α2h3L,αγ1/2h3/2N1/2+α2h3LN2,α2h2N} | (4.28) |
which concludes the proof.
Remark 4.5. We remark that the first inequality in (4.26) corresponds to the choice ℓ:=αh3/2/(γ1/2n3/2). Therefore in the regime in which the energy scales as αγ1/2h3/2N1/2, delamination necessarily occurs over a length which is (up to a factor) equal to ℓ:=αh3/2/(γ1/2n3/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
h3L[|R0+−R−α|2+|R0−−Rα|2]≤cEh[u]. | (4.29) |
We recall that R0+−R−α and R0−−Rα are the differences in rotation inside the two non-delaminated parts of the sample. Therefore if the energy is significantly smaller than cα2h3/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 γ 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 γ.
The bending moment is displayed in Figure 6. We expect that the moment is continuous in the bending angle α, 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 γ). 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
![]() |
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 |