LDG approximation of large deformations of prestrained plates

https://doi.org/10.1016/j.jcp.2021.110719Get rights and content

Highlights

  • Formal derivation of prestrained hyper-elastic plates model.

  • Novel reformulation amenable to variational approximation methods.

  • Novel finite element methods based on reconstructed Hessians.

  • Insightful numerical experiments.

Abstract

A reduced model for large deformations of prestrained plates consists of minimizing a second order bending energy subject to a nonconvex metric constraint. The former involves the second fundamental form of the middle plate and the latter is a restriction on its first fundamental form. We discuss a formal derivation of this reduced model along with an equivalent formulation that makes it amenable computationally. We propose a local discontinuous Galerkin (LDG) finite element approach that hinges on the notion of reconstructed Hessian. We design discrete gradient flows to minimize the ensuing nonconvex problem and to find a suitable initial deformation. We present several insightful numerical experiments, some of practical interest, and assess various computational aspects of the approximation process.

Introduction

Natural and manufactured phenomena abound where thin materials develop internal stresses, deform out of plane and exhibit nontrivial 3d shapes. Nematic glasses [30], [31], natural growth of soft tissues [22], [38] and manufactured polymer gels [25], [26], [37] are chief examples. Such incompatible prestrained materials may be key constituents of micro-mechanical devices and be subject to actuation. A model postulates that these plates may reduce internal stresses by undergoing large out of plane deformations u as a means to minimize an elastic energy E[u] that measures the discrepancy between a reference (or target) metric G and the orientation preserving realization u of it. The strain tensor ϵG(u), given byϵG(u):=12(uTuG), measures such discrepancy and yields the following elastic energy functional for prestrained isotropic materials in a 3d reference body B and without external forcingE[u]:=Bμ|G12ϵG(u)G12|2+λ2tr(G12ϵG(u)G12)2, where μ,λ are the Lamé constants [17], [18], [36]. A deformation u:BR3 such that ϵG(u)=0 is called isometric immersion. If such a map exists, then the material can attain a stress-free equilibrium configuration, i.e., E[u]=0. However, the existence of an isometric immersion u of class H2(B) for any given smooth metric G is not guaranteed in general and it constitutes an outstanding problem in differential geometry. In the absence of such a map, the infimum of E[u] is strictly positive and the material has a residual stress at free equilibria.

Slender elastic bodies are of special interest in many applications and our main focus. In this case, the 3d domain B can be viewed as a tensor product of a 2d domain Ω, the midplane, and an interval of length s, namely Ω×(s2,s2). Developing dimensionally-reduced models as s0 is a classical endeavor in nonlinear elasticity. Upon rescaling E[u] with a factor of the form sβ, several 2d models can be derived in the limit s0. A geometrically nonlinear reduced energy was obtained formally by Kirchhoff in his seminal work of 1850. An ansatz-free rigorous derivation for isotropic materials was carried out in the influential work of Friesecke, James and Müller in 2002 [19] via Γ-convergence for β=3. This corresponds to the bending regime of the nonlinear Kirchhoff plate theory.

If the target metric G is the identity matrix, there is no in-plane stretching and shearing of the material leaving bending as the chief mechanism of deformation; an excellent example examined in [19] is the bending of a sheet of paper. For a generic metric G that does not depend on s and is uniform across the thickness, Efrati, Sharon and Kupferman derived a 2d energy which decomposes into stretching and bending components [17]; the former scales linearly in s whereas the latter does it cubically. The first fundamental form of the midplane characterizes stretching while the second fundamental form accounts for bending. The thickness parameter s appears in the reduced energy and determines the relative weight between stretching and bending.

The asymptotic limit s0 requires a choice of scaling exponent β. The bending regime β=3 has been studied by Lewicka and collaborators [29], [8], while [20], [7], [27], [28] discussed other exponents β. For instance, β=5 corresponds to the Föppl von Kárman plate theory, which is suitable for moderate deformations. Different energy scalings select specific asymptotic relations between the prestrain metric G and deformations u. For instance, for β=3 and metrics G of the formG(x,x3)=G(x)=[g(x)001]xΩ,x3(s/2,s/2), with g:ΩR2×2 a symmetric uniformly positive definite matrix, considered for instance in [29], [17], the first fundamental form

of parametrizations y:ΩR3 of the midplane must satisfy the following pointwise metric constraint as s0 this accounts for the stretching and shearing of the midplane. Moreover, the scaled elastic energy s3E[u] turns out to Γ-converge to the reduced bending energy which depends solely on the second fundamental form
of y in the absence of external forcing [29], [8]. It is known that E[y]>0 provided that the Gaussian curvature of the surface y(Ω) does not vanish identically [29], [8]. We illustrate this in Fig. 1.

In this article, we present a numerical study of the minimization of (5) subject to the constraint (4) with either Dirichlet or free boundary conditions. We start in Section 2 with a justification of (2) followed by a formal derivation of (4) and (5) as the asymptotic limit of s3E[u] as s0 when the 3d target metric is of the form (3). Moreover, we show an equivalent formulation that basically replaces the second fundamental form

by the Hessian D2y, which makes the constrained minimization problem amenable to computation. This derivation is, however, trickier than that in [3], [12] for single layer plates and [5], [4], [11] for bilayer plates.

The numerical treatment of the ensuing fourth order problem is a challenging and exciting endeavor. In [3], [5], [4], the discretization hinges on Kirchhoff elements for isometries y, i.e., g=I2R2×2 is the identity matrix. The approximation of y in [12], [11] relies on a discontinuous Galerkin (dG) method. In all cases, the minimization problem associated with the nonconvex constraint (4) resorts to a discrete H2-gradient flow approach. In addition, a Γ-convergence theory is developed in [3], [5], [12]. In this paper, inspired by [15], we design a local discontinuous Galerkin method (LDG) for gI2 that replaces D2y by a reconstructed Hessian Hh[yh] of the discontinuous piecewise polynomial approximation yh of y. Such discrete Hessian Hh[yh] consists of three distinct parts: the broken Hessian Dh2yh, the lifting Rh([hyh]) of the jump of the broken gradient hyh of yh, and the lifting Bh([yh]) of the jumps of yh itself. Lifting operators were introduced in [6] and analyzed in [13], [14]. The definition of Rh and Bh is motivated by the liftings of [32], [33] leading to discrete gradient operators.

It is worth pointing out prior uses of Hh[yh]. Discrete Hessians were instrumental to study convergence of dG for the bi-Laplacian in [35] and plates with isometry constraint in [12]. In the present contribution, Hh[yh] makes its debut as a chief constituent of the numerical method. We introduce Hh[yh] in Section 3 along with the LDG approximation of (5) and the metric defect Dh[yh] that relaxes (4) and makes it computable. We also discuss two discrete H2-gradient flows, one to reduce the bending energy (5) starting from yh0, and the other to diminish the stretching energy and make Dh[yh0] as small as possible. The former leads to Algorithm 1 (gradient flow) and the latter to Algorithm 2 (initialization) of Section 3. We reserve Section 4 for implementation aspects of Algorithm 1, Algorithm 2. We present several numerical experiments, some of practical interest, in Section 5 that document the performance of the LDG approach and illustrate the rich variety of shapes achievable with the reduced model (4)-(5). We close the paper with concluding remarks in Section 6.

Section snippets

Problem statement

Let Ωs:=Ω×(s/2,s/2)R3 be a three-dimensional plate at rest, where s>0 denotes the thickness and ΩR2 is the midplane. Given a Riemannian metric G:ΩsR3×3 (symmetric uniformly positive definite matrix), we consider 3d deformations u:ΩsR3 driven by the strain tensor ϵG(u) of (1) that measures the discrepancy between uTu and G; hence, the 3d elastic energy E[u]=0 whenever ϵG(u)=0. We say that G is the reference (prestrained or target) metric. An orientable deformation u:ΩsR3 of class H2(Ω)

Numerical scheme

We propose here a local discontinuous Galerkin (LDG) method to approximate the solution of the problem (26). LDG is inspired by, and in fact improves upon, the previous dG methods [11], [12] but they are conceptually different. LDG hinges on the explicit computation of a discrete Hessian Hh[yh] for the discontinuous piecewise polynomial approximation yh of y, which allows for a direct discretization of Eh[yh] in (36), including the trace term. We refer to the companion paper [9] for a

Implementation

We make a few comments on the implementation of the gradient flow (61)-(62), built in Algorithm 1, and the resulting linear algebra solver used at each step.

Numerical experiments

In this section, we present a collection of numerical experiments to illustrate the performance of the proposed methodology. We consider several prestrain tensors g, as well as both ΓD (Dirichlet boundary condition) and ΓD= (free boundary condition). Algorithm 1, Algorithm 2 are implemented using the deal.ii library [2] and the visualization is performed with paraview [1]. The color code is the following: (multicolor figures) dark blue indicates the lowest value of the deformation's third

Conclusions

In this article, we design and implement a numerical scheme for the simulation of large deformations of prestrained plates. Our contributions are:

1. Model and asymptotics. We present a formal asymptotic limit of a 3d hyperelastic energy in the bending regime. The reduced model, rigorously derived in [8], consists of minimizing a nonlinear energy involving the second fundamental form of the deformed plate and the target metric under a nonconvex metric constraint. We show that this energy is

CRediT authorship contribution statement

All the authors contributed to all the aspects of this work.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Ricardo H. Nochetto and Shuo Yang were partially supported by the NSF Grants DMS-1411808 and DMS-1908267.

Andrea Bonito and Diane Guignard were partially supported by the NSF Grant DMS-1817691.

References (38)

  • K. Bhattacharya et al.

    Plates with incompatible prestrain

    Arch. Ration. Mech. Anal.

    (2016)
  • A. Bonito et al.

    Numerical analysis of the LDG method for large deformations of prestrained plates

  • A. Bonito et al.

    Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method

    SIAM J. Numer. Anal.

    (2010)
  • A. Bonito et al.

    DG approach to large bending deformations with isometry constraint

    Math. Models Methods Appl. Sci.

    (2021)
  • F. Brezzi et al.

    Discontinuous finite elements for diffusion problems

  • F. Brezzi et al.

    Discontinuous Galerkin approximations for elliptic problems

    Numer. Methods Partial Differ. Equ.

    (2000)
  • B. Cockburn et al.

    The local discontinuous Galerkin method for time-dependent convection-diffusion systems

    SIAM J. Numer. Anal.

    (1988)
  • M.P. do Carmo

    Differential Geometry of Curves and Surfaces

    (1976)
  • E. Efrati et al.

    Hyperbolic non-euclidean elastic strips and almost minimal surfaces

    Phys. Rev. E

    (2011)
  • Cited by (14)

    • An anisotropic adaptive method for the numerical approximation of orthogonal maps

      2022, Journal of Computational and Applied Mathematics
      Citation Excerpt :

      The solution being piecewise linear, the key aspect is to approximate the discontinuity lines (the so-called singular set), and this has been the reason to introduce anisotropic adaptive methods. Besides the original origami application [4], such orthogonal maps formulation is or will be used for some applications, e.g., deformations of prestrained plates [14,15], aerospace engineering [16], cartography [17], differential geometry [18] and, to some extend, robotics [19,20]. More recently, it has been applied to several fields of material science, such as 3D printing and material composites [21–24] (e.g. for self folding materials).

    • Error estimates for a linear folding model

      2024, IMA Journal of Numerical Analysis
    • A nonlinear bending theory for nematic LCE plates

      2023, Mathematical Models and Methods in Applied Sciences
    • A Homogenized Bending Theory for Prestrained Plates

      2023, Journal of Nonlinear Science
    View all citing articles on Scopus
    View full text