Formation of wrinkles on a coated substrate using manifold‐valued finite elements

This article treats finite element simulations of controlled wrinkle formation experiments of a soft bulk material with a thin, stiff layer on top. The wrinkling process is triggered by a stress mismatch between the bulk material and the thin layer. For the finite element simulations, we model the bulk material using a three‐dimensional hyperelastic material and the thin layer with a geometrically nonlinear elastic Cosserat shell. For the finite element simulations, we model the bulk material using a three‐dimensional hyperelastic material and the thin layer with a geometrically nonlinear elastic Cosserat shell. We use Lagrange finite elements for the bulk material and geodesic finite elements for the shell. The resulting minimization problem is nonlinear and nonconvex. We prove existence of minimizers in the continuous and the discrete function space. Finally, we solve the resulting nonconvex minimization problem numerically using a Riemannian trust‐region algorithm and compare our simulations to real experiments.


INTRODUCTION
This article treats finite element simulations of controlled wrinkle formation experiments on a coated substrate conducted at the Leibnitz Institute for Polymer Research by Knapp et al. [1].The systems used for the experiments consist of a soft polydimethylsiloxane (PDMS) layer with a thin, stiff layer on top.The wrinkling process is triggered by a stress mismatch between the bulk and the thin layer.To create the stress mismatch, the bulk material is first uni-axially stretched and then the thin layer is created by a low-pressure plasma treatment in a vacuum chamber, with the gases  2 or  2 .Under subsequent relaxation, wrinkles form as shown in Figure 1.Their wavelength depends on the choice of the process gas and the plasma process parameters such as the duration of the treatment.The use of thin silicon masks placed directly on the PDMS allows to sharply restrict the plasma-exposed area.Sequential exposures of the same sample to different treatment processes with and without a mask allow to locally modify the layer thickness and stiffness.This allows to locally control the wavelength of the resulting wrinkles and to trigger the formation of branches and line defects located at the boundary between areas of different wavelengths.
While closed-form expressions of the wrinkle wavelength are available for some simple cases, scenarios including branching and line defects require numerical simulations.
F I G U R E 1 Schematic illustration of the masked wrinkling process [1].

F I G U R E 2
Sketch of the mathematical model.The stress-free configuration of the shell is determined by the initial stretch of the bulk.

MATHEMATICAL MODEL
To model the system, we assume that the deformation process is time-independent and elastic.We use a Mooney-Rivlin material model for the bulk material and the Cosserat shell model from [2,3] for the thin layer.Consider an open and bounded Lipschitz domain Ω ⊂ ℝ 3 describing the undeformed, stress-free state of the bulk material.Let the boundary Ω be split into a Dirichlet part Γ  ⊂ Ω and a Neumann part Γ  ⊂ Ω, where Γ  ∩ Γ  = ∅, and Γ  ∪ Γ  = Ω.Futhermore, let a part Γ  ⊂ Γ  be diffeomorphic to a domain  ℎ ⊂ ℝ 2 , for a diffeomorphism  0 ∈  ∞ ( ℎ , Γ  ).
For the wrinkling process, we first consider Ω to be deformed.Denote the deformation function by  0 , the deformed state by and assume that ( 0 ) | Γ  ∈  2 (Γ  , ℝ 3 ).Then, we consider the deformed bulk material to be coupled with the Cosserat shell as shown in Figure 2.For the Cosserat shell, denote the stress-free state of the shell surface by   , where   is a two-dimensional manifold given by the parametrization The stress-free state of the Cosserat shell surface   and the stress-free state Ω of the bulk material do not match.This stress mismatch results in the wrinkle formation  w ∶ Ω 0 → ℝ 3 , as shown in Figure 2. Denote the deformation function mapping Ω to  w (Ω 0 ) by  and the deformation function for the Cosserat shell mapping  ℎ to (Γ  ) by .The deformation of the shell is then for 2 ) ∈  ℎ .Cosserat models introduce an orientation of the material particles in addition to their position, so the current state of the shell surface is not only described by the deformation function  but also by an independent microrotation function denoted by  ∶  ℎ → SO (3).
The energy functional  for the coupled system is It sums the energy density function  bulk (∇) of the Mooney-Rivlin material model of second order from [4] and the energy density function of the Cosserat shell  shell (  ,   ) =  memb (  ,   ) +  bend (  ) from [2,3].The definitions of the shell strain tensor   and the shell bending-curvature tensor   can be found in [2,3].We want to find a minimizer of the functional (, ) in the admissible set where  * ∈  1, (Ω, ℝ 3 ) is a given Dirichlet function with  * | Γ  ∈  1 (Γ  , ℝ 3 ) and the manifold-valued Sobolev space  1 (Γ  , SO(3)) is defined as in [5].
Theorem 1. Assume that: (A1) The parametrization  0 = ( 0 ) | Γ  • 0 of the stress-free state of the Cosserat shell surface fulfills Then there is a pair of a deformation and a microrotation function in the admissible set  that minimizes the functional (2).

FINITE ELEMENT SIMULATION
Assume that there is a triangulation  3 of the domain Ω which induces a matching triangulation  2 of  ℎ .The codomain space of the deformation function  is the linear vector space ℝ 3 .Therefore, the deformation can be approximated with Lagrange finite element functions.Denote the Lagrange finite element space of order   by  poly   (Ω, ℝ 3 ).To discretize the microrotation function  introduced by the Cosserat shell model, Lagrange finite element functions cannot be used because microrotations map into the nonlinear Riemannian manifold SO(3).Instead, we use geodesic finite elements [3,5].Denote the geodesic finite element space of order   by  geo   ( ℎ , SO(3)).As shown experimentally in [3] and [5], no locking occurs for this type of shell discretization if finite elements of second order are chosen.

Theorem 2. Assume that:
(A1) The parametrization  0 of the stress-free state of the Cosserat shell surface fulfills

} contains at least one function pair with finite energy. (A3) -(A5) From Theorem 1 hold.
Then there is a deformation and microrotation function pair in the admissible set  ℎ that minimizes the functional (2).
Proof.[7,Section 3.3] □ The minimization problem on the set  ℎ then yields an algebraic minimization problem which we solve using a Riemannian trust-region method from [8] and a monotone multigrid method from [9,10] for the inner problems.
We perform a simulation of an experiment from [1].The bulk is a rectangular block of 60 × 20 × 20 μm 3 .The stiff layer covers a region measuring 30 × 20 μm 2 as shown in Figure 3.This block is discretized by a cubical grid which is manually refined in the vicinity of the shell yielding a grid as shown in [3].In total, the grid contains 15 642 vertices and 79 578 second-order Lagrange nodes.The shell part contains 3 993 vertices and 9 801 second-order Lagrange nodes.The parameters   of the Mooney-Rivlin material from [4] were derived from uniaxial tensile tests with the material used in the experiment (Sylgard 184 with a base to hardener component ratio of 2:1).The material parameters are  01 = 1.94 MPa,  10 = −1.67MPa,  02 = 6.52 MPa,  20 = 2.42 MPa,  11 = −7.34MPa, and  = 57 MPa.
In the experiment the base material was treated with  2 plasma for 60  and 240 s as shown in Figure 3.The thickness and Young's modulus  of the differently treated areas were measured by a quantitative nanoscale mechanical characterization and we assumed a Poisson's ratio of  = 0.335.From that, we calculated the Lamé parameters  and  for the layer.All material parameters for the shell are listed in Table 1.
The simulation shown in Figure 4 was implemented in C++ based on the Dune libraries [11].It was executed using 24 tasks of a Intel(R) Xeon(R) CPU E5-2680 with 2.50 GHz and 64 GB RAM.It took 56 trust-region steps and 24.35 h to solve the problem.
We compare the experiment and the simulation in terms of the wavelength, amplitude and shape of the wrinkles in the differently treated regions.For this, we investigate a square area of 18 × 18 μm 2 of the treated sample by an atomic force microscope (AFM) scan and post-process the resulting data, yielding the top right photograph of Figure 4. To compare our simulation with the experiment, we analogously investigate a square area of 18 × 18 μm 2 of the simulation result, yielding the top left image of Figure 4.
Overall, the wavelengths as well as the amplitudes measured from the experiment and from the simulation and calculated match well, as shown in Table 2.For both the simulation and the experiment, branching occurs in the transition zone between the adjacent wrinkling areas, as shown in Figure 4.For a more thorough comparison, we refer to [1].
A2) The admissible set  from Equation (3) contains at least one function pair with finite energy.(A3) The bulk energy density function  bulk is polyconvex and bounded from below, in the sense that there are constants  > 0 and  > −∞ such that  bulk (∇) ≥ |∇|  + , for all  ∈  and  > 3. (A4) The energy density function  bulk fulfills  bulk (∇) → ∞ if det(∇) → 0. (A5) The coefficients of the energy density functions  memb and  bend of the shell fulfill the assumptions of[6, Theorem 3.3].

F I G U R E 3
The undeformed grid for the bulk.The top view shows the two areas where plasma treatment has been applied for 240 and 60 s, respectively, to create surface layers of different thickness and stiffness.TA B L E 1 Material parameters for the Cosserat shell.

F I G U R E 4
Comparison of the simulation and the experimental wrinkle scenario: In the two top pictures we see a square area of 18 μm × 18 μm from the center of Γ  .In the two bottom pictures we see the comparison of the branching area[1].TA B L E 2Measured thickness and stiffness values for the differently treated areas and resulting wavelengths and amplitudes.