Modeling and simulation of thermally actuated bilayer plates

TL;DR

This paper develops a mathematical model for thermally actuated polymer bilayer plates, enabling simulation of large bending deformations due to thermal stimuli, with practical discretization and iterative solution methods.

Contribution

It introduces a new nonlinear bending model for thermally responsive bilayer plates and provides a practical simulation scheme for folding behaviors.

Findings

Successfully simulates large thermal bending deformations.

Demonstrates folding of practical compliant structures.

Provides a discretization method suitable for engineering applications.

Abstract

We present a mathematical model of polymer bilayers that undergo large bending deformations when actuated by non-mechanical stimuli such as thermal effects. The simple model captures a large class of nonlinear bending effects and can be discretized with standard plate elements. We devise a fully practical iterative scheme and apply it to the simulation of folding of several practically useful compliant structures comprising of thin elastic layers.

Equations

$$W(E):=\mu |E{|}^2+\frac{\lambda }{2}(\mathrm{tr}E{)}^2,$$

$$W(E_{\Theta })=W(E)-(2\mu +3\lambda )\alpha \Theta \mathrm{tr}E+c_{\Theta }$$

$$\mu|E_{\Theta}|^{2}\leq W(E_{\Theta})\leq\big{(}\mu+\frac{3\lambda}{2}\big{)}|E_{\Theta}|^{2}.$$

$$W_{\Theta}(F):=\frac{\mu}{4}|E_{\Theta}|^{2}=\frac{\mu}{4}\big{|}F^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}F-(1+\alpha\Theta)I\big{|}^{2}.$$

$$W_{\Theta}(\nabla\mathbf{u},\mathbf{x},t)=\frac{\mu(\mathbf{x})}{4}\Big{|}\nabla\mathbf{u}^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla\mathbf{u}-(1+\alpha(\mathbf{x})\Theta(\mathbf{x},t))I\Big{|}^{2},$$

$$\alpha (\mathbf{x}):=\{\begin{array}{ll}+\alpha & \text{mbox}for{x}_3>0,\\ -\alpha & \text{mbox}for{x}_3<0,\end{array}\mu (\mathbf{x}):=\{\begin{array}{ll}{\mu }^{+}(x^{\prime })& \text{mbox}for{x}_3>0,\\ {\mu }^{-}(x^{\prime })& \text{mbox}for{x}_3<0.\end{array}$$

$$\mathbf{u}({\mathbf{x}}^{\prime },x_3,t)=\mathbf{y}({\mathbf{x}}^{\prime },t)+x_3\mathbf{b}({\mathbf{x}}^{\prime },t)$$

$$\nabla \mathbf{u}(\mathbf{x})=[{\nabla }^{\prime }\mathbf{y}({\mathbf{x}}^{\prime }),\mathbf{b}({\mathbf{x}}^{\prime })]+x_3[{\nabla }^{\prime }\mathbf{b}({\mathbf{x}}^{\prime }),0]$$

$$I_{\Theta }[\mathbf{u}]=\frac{1}{{\delta }^3}{\int }_{{\omega }_{\delta }}{W}_{\Theta }(\nabla \mathbf{u},\cdot )$$

$$\begin{split}I_{\Theta}[\mathbf{u}]=\frac{1}{\delta^{3}}\int_{\omega}\int_{-\delta/2}^{\delta/2}\frac{\mu}{4}\Bigg{\{}\Bigg{|}&\begin{bmatrix}(\nabla^{\prime}\mathbf{y})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla\mathbf{y}-(1\pm\alpha\Theta)I_{2}&0\\ 0&|\mathbf{b}|^{2}-(1\pm\alpha\Theta)\end{bmatrix}\\ +x_{3}&\begin{bmatrix}(\nabla^{\prime}\mathbf{b})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{y}+(\nabla^{\prime}\mathbf{y})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{b}&(\nabla^{\prime}\mathbf{b})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\mathbf{b}\\ \mathbf{b}^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{b}&0\end{bmatrix}\\ +x_{3}^{2}&\begin{bmatrix}(\nabla^{\prime}\mathbf{b})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{b}&0\\ 0&0\end{bmatrix}\Bigg{|}^{2}\Bigg{\}}.\end{split}$$

$$G=(\nabla^{\prime}\mathbf{y})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{y},\quad H=-(\nabla^{\prime}{\boldsymbol{\nu}})^{\rotatebox[origin={c}]{180.0}{\footnotesize$\perp$}}\nabla^{\prime}\mathbf{y}.$$

$$[{\nabla }^{\prime }\mathbf{y}{]}^T{\nabla }^{\prime }\mathbf{y}=I_2,$$

$$\delta^{-3}\int_{-\delta/2}^{\delta/2}\frac{\mu}{4}|2x_{3}H\mp\alpha\Theta I_{2}|^{2}=\frac{\overline{\mu}}{12}\Big{(}|H-\overline{\alpha}\theta I_{2}|^{2}-\frac{4}{3}\overline{\alpha}^{2}\theta^{2}\Big{)},$$

$$\overline{\alpha }:=\underset{\delta \to 0}{\lim }\frac{3\alpha }{\delta }.$$

$$E[\mathbf{y}]=\frac{1}{12}\int_{\omega}\overline{\mu}\big{|}H-\overline{\alpha}\theta I_{2}\big{|}^{2},$$

$$\overline{\sigma }:=\frac{1}{2}({\sigma }^{+}+{\sigma }^{-}),\overline{\kappa }:=\frac{1}{2}({\kappa }^{+}+{\kappa }^{-})$$

$$\nabla \Theta ={\nabla }_{\gamma }\Theta +{\partial }_{\mathbf{\nu }}\Theta \mathbf{\nu },$$

$$\displaystyle\nabla_{\gamma}\Theta(\mathbf{x}^{\prime},x_{3})=\big{(}I-{\boldsymbol{\nu}}(\mathbf{x}^{\prime})\otimes{\boldsymbol{\nu}}(\mathbf{x}^{\prime})\big{)}\nabla\Theta(\mathbf{x}^{\prime},x_{3}),$$

$${\partial }_{\mathbf{\nu }}\Theta ({\mathbf{x}}^{\prime },x_3)=\mathbf{\nu }({\mathbf{x}}^{\prime })\cdot \nabla \Theta ({\mathbf{x}}^{\prime },x_3),$$

$${\partial }_{\mathbf{\nu }}\Theta ({\mathbf{x}}^{\prime },x_3)\to 0\text{as }\delta \to 0.$$

$$\frac{d}{dt}{\int }_{{\gamma }_{\delta }(t)}\sigma \Theta =-{\int }_{\partial {\gamma }_{\delta }(t)}\mathbf{q}\cdot \mathbf{n}=-{\int }_{{\gamma }_{\delta }(t)}\mathrm{div}\mathbf{q}.$$

$$\frac{d}{dt}\int_{\gamma_{\delta}(t)}\sigma\Theta=\int_{\gamma_{\delta}(t)}\sigma\big{(}\partial_{t}\Theta+\nabla\Theta\cdot\partial_{t}\mathbf{u}+\Theta\operatorname{div}\partial_{t}\mathbf{u}\big{)}.$$

$$\frac{1}{\delta }{\int }_{{\gamma }_{\delta }(t)}\sigma \nabla \Theta \cdot {\partial }_t\mathbf{u}=\frac{1}{\delta }{\int }_{{\gamma }_{\delta }(t)}\sigma {\nabla }_{\gamma }\Theta \cdot {\partial }_t\mathbf{u}+\sigma {\partial }_{\mathbf{\nu }}\Theta \mathbf{\nu }\cdot {\partial }_t\mathbf{u}\to {\int }_{\gamma (t)}\overline{\sigma }{\nabla }_{\gamma }\theta \cdot {\partial }_t\mathbf{y},$$

$${\partial }_{x_3}\mathbf{F}(\mathbf{x})=\nabla \mathbf{f}(\mathbf{u}(\mathbf{x})){\partial }_{x_3}\mathbf{u}(\mathbf{x})=\nabla \mathbf{f}(\mathbf{u}(\mathbf{x}))\mathbf{b}({\mathbf{x}}^{\prime })\forall \mathbf{x}\in {\gamma }_{\delta }.$$

$$\mathrm{div}\mathbf{f}={\mathrm{div}}_{\gamma }\mathbf{f}+\mathbf{\nu }\cdot \nabla \mathbf{f}\mathbf{\nu }={\mathrm{div}}_{\gamma }\mathbf{f}+|\mathbf{b}{|}^{-2}\mathbf{b}\cdot {\partial }_{x_3}\mathbf{F}.$$

$$\mathrm{div}{\partial }_t\mathbf{u}={\mathrm{div}}_{\gamma }{\partial }_t\mathbf{u}.$$

$$\frac{d}{dt}\frac{1}{\delta}\int_{\gamma_{\delta}(t)}\sigma\Theta\to\int_{\gamma}\overline{\sigma}\big{(}\partial_{t}\theta+\nabla_{\gamma}\theta\cdot\partial_{t}\mathbf{y}+\theta\operatorname{div}_{\gamma}\partial_{t}\mathbf{y}\big{)}=\frac{d}{dt}\int_{\gamma(t)}\overline{\sigma}\theta.$$

$$\mathrm{div}\mathbf{q}={\mathrm{div}}_{\gamma }\mathbf{q}+{\mathrm{div}}_{\gamma }^{\perp }\mathbf{q},$$

$${\int }_{{\gamma }_{\delta }}{\mathrm{div}}_{\gamma }^{\perp }\mathbf{q}=0.$$

$$\frac{1}{\delta }{\int }_{{\gamma }_{\delta }}{\mathrm{div}}_{\gamma }\kappa \nabla \Theta =\frac{1}{\delta }{\int }_{{\gamma }_{\delta }}{\mathrm{div}}_{\gamma }\kappa \nabla \gamma \Theta +\frac{1}{\delta }{\int }_{{\gamma }_{\delta }}{\mathrm{div}}_{\gamma }\kappa ({\partial }_{\mathbf{\nu }}\Theta \mathbf{\nu })\to {\int }_{\gamma (t)}{\mathrm{div}}_{\gamma }\overline{\kappa }{\nabla }_{\gamma }\theta ,$$

$${\mathrm{div}}_{\gamma }({\partial }_{\mathbf{\nu }}\Theta \mathbf{\nu })={\nabla }_{\gamma }({\partial }_{\mathbf{\nu }}\Theta )\cdot \mathbf{\nu }+{\partial }_{\mathbf{\nu }}\Theta {\mathrm{div}}_{\gamma }\mathbf{\nu }={\partial }_{\mathbf{\nu }}\Theta \mathrm{tr}H\to 0\text{as }\delta \to 0.$$

Full text

A.B. was partially supported by NSF Grant DMS-1254618 and AFOSR Grant FA9550-14-1-0234.

A.H.M. was partially supported by AFOSR Grant FA9550-14-1-0234.

R.H.N was partially supported by NSF Grant DMS-1411808, Institut Henri Poincaré and the Simons Visiting Professorship (Oberwolfach)

1. Introduction

Bilayer polymers are appealing for the development of autonomous lightweight foldable structures, such as drug delivery vesicles, flexible constructions for soft robots, self-deployable sun sails in spacecraft, and morphing structures, since they can be manufactured into various shapes with tunable material properties for each layer. A characteristic feature is that they can undergo controlled large deformations via relatively small external stimuli. Bilayers having two polymeric layers with different thermo-responsive and swelling characteristics have been shown capable of forming various practically useful folding shapes [28, 6, 2]. Different expansion/contraction characteristics of the two layers when exposed to temperature changes or solvent diffusions lead to out of plane rotations and curvature changes of the bilayers which trigger folding, cf. Fig. 1. When the two layers have significantly different expansion/contraction behaviors, often characterized by the coefficient of thermal or moisture expansion, folding can be achieved with relatively smaller stimuli. However, high contrast in mechanical properties, i.e., elastic modulus, between the two layers can induce high stress discontinuities between the layers, leading to delamination. Another mechanism for folding and bending of bilayers can be achieved by integrating two electro-active polymers with the opposite directions of through-thickness poling axes and when the bilayers are subjected to electric potential through the thickness one layer would expand while the other would contract [31, 30, 21]. One of the advantages of combining two polymers with different responsive characteristics with regard to their non-mechanical performances is that the mechanical properties of polymers do not vary significantly, e.g., extensional elastic moduli of various polymers are typically between $0.5-5.0\,{\rm GPa}$, which can minimize stress discontinuities at the interfaces between the two layers and thus can avoid failure due to delamination. Alternatively, bilayers can be formed by combining two different types of materials, i.e., metal- or ceramic-polymer, which have significant differences in their mechanical and non-mechanical properties [25, 14].

We are specifically interested in bilayers comprising of two polymeric layers that can undergo large deformations when exposed to non-mechanical stimuli, such as temperature changes or fluid sorption. Due to the slender nature of the bilayers, large deformations are mainly governed by rotations while the strains in the bilayer are relatively small [26, 30], thereby leading to negligible stretching and transverse shear effects. Furthermore, the mid-surface of the bilayer plate is considered to be inextensible and non-shearable, thereby leaving bending as a chief mechanism for shape deformation. These basic mechanical assumptions lead to a reduced Kirchhoff plate model for the deformation with a preferred curvature tensor and an isometry constraint. The former encodes the mismatch between bilayers while the latter reflects the property that distances among points on the mid surface do not change with shape deformation. In addition, we ignore inertial effects and assume a quasi-static evolution of the plate.

The bilayer plate is driven by heat conduction or fluid sorption. We model this diffusion process with a linear Fourier heat conduction or Fickian diffusion law. Since the deformation is an isometry, the diffusion equation is insensitive to the plate deformation, whereas the temperature affects the mismatch between lower and upper layers of the plate. In several experimental studies, the bilayers are insulated on top and bottom and subjected to uniform environmental conditions on their sides, e.g., at dry conditions the bilayer is immersed in fluid or from curing at high temperature the bilayer is cooled down to room temperature or vice versa. In such situations, the entire polymers will be at uniform temperature or will have uniform fluid content at the steady state and the diffusion process occurs from the boundaries of the bilayers. We model this with homogeneous Neumann boundary condition on top and bottom of the plate, and either Neumann or Robin boundary conditions on the sides of the plate. We derive a reduce diffusion model that accounts for these effects.

Our investigations are motivated by lab experiments reporting problems in the controlled fabrication of nanotubes, cf. [28, 35, 33, 34]. In particular, rectangular bilayer plates occasionally start bending from the corners leading to formation of so-called dog-ears, thus failing to attain the desired cylindrical shape. Our mathematical model and corresponding simulations indicate that this may be attributed to too rapid changes in the environment that yield high concentrations of the diffusing quantity at the corners of the bilayer. Moreover, we present computational experiments resembling origami structures such as self-assembling cube [29, 13] and deployable airfoil [17], as well as particle encapsulation which is of interest in drug targeting [27]. Although additional effects become relevant at the nanoscale, see e.g. [32], our computational model could be utilized as a simple tool to explore and predict conditions and configurations that enable a controlled production of nanoscale and microscale devices. The method is advantageous with respect to direct simulation of full three-dimensional thermo-elasticity in terms of simplicity of implementation, efficiency, performance, and flexibility.

Our simple two-dimensional model is thus able to capture the principal mechanical and thermal effects responsible for large deformations observed in lab experiments of thermally actuated bilayer plates. This suffices for the aforementioned applications. In fact, the present model and analyses can help designers in simulating desired shape changes and determining external stimuli to be prescribed prior to fabricating flexible bilayer systems. We refer the reader to [1, 6, 15, 18, 23, 24, 25, 27, 28, 36] for further discussion on the applications. In contrast, a three-dimensional model capturing strong swelling effects and which leads to a two-way coupling of the diffusion and deformation processes has been investigated in [20].

Most of available studies with regards to folding of polymeric composite structures have been on fabrication and experiments. The mathematical description and numerical treatment of large bending deformations of elastic solids has recently undergone some important development. Dimensionally reduced models have been derived rigorously from three-dimensional hyperelasticity [10, 21] and numerical methods capable of approximating large rotations correctly have been devised and analyzed in [3, 4, 5]. In this work we extend the approach from [5] by including temperature dependence in the preferred curvature tensor. The mechanical part consists of a nonlinear Kirchhoff model for which we devise a finite element discretization based on standard plate elements. The heat equation decouples from the mechanical equation due to the inextensibility of the plate and we approximate it with standard finite element methods together with the backward Euler time stepping algorithm. The proposed iterative numerical method for the coupled system is roughly ten times faster than the one used in [5].

The outline of this article is as follows. In Section 2 we introduce the thermo-mechanical mathematical model and describe its dimension reduction. Corresponding partial differential equations are formulated in Section 3. We then present in Section 4 the temporal and spatial discretization of the nonlinear and constrained system composed of second order diffusion and fourth order bending equations. We report several intruiguing and practically useful numerical experiments in Section 5.

2. Mathematical model

2.1. Hyperelastic materials

We model polymers as isotropic and elastic materials. According to St. Venant–Kirchhoff description, the mechanical behavior of the system is governed by a hyperelastic stored energy density

[TABLE]

where E:=F^{\rotatebox[origin={c}]{180.0}{\footnotesize\perp}}F-I is the Green–Lagrange strain tensor related to the deformation gradient $F$, and $\lambda,\mu$ are the (first and second) Lamé constants [19]; hereafter we let $|E|^{2}=E:E$ be the square of the Frobenius norm of $E$ and ${\rm tr}E=E:I$ be the trace of $E$. As strain tensor we use the temperature dependent quantity E_{\Theta}:=F^{\rotatebox[origin={c}]{180.0}{\footnotesize\perp}}F-(1+\alpha\Theta)I\in\mathbb{R}^{3\times 3} which generalizes the Green–Lagrange strain tensor. Here $\alpha,\Theta$ are the thermal expansion coefficient and temperature change. Therefore, a simple calculation yields

[TABLE]

with $c_{\Theta}:=(3\mu+\frac{9}{2}\lambda)\alpha^{2}\Theta^{2}$. We note that the elastic material constants $\lambda,\mu$ can also change with temperature or fluid concentration. We simplify the model further upon realizing that $({\rm tr}E)^{2}\leq 3|E|^{2}$, whence

[TABLE]

This equivalence motivates our choice of energy

[TABLE]

which in turn allows us to follow the arguments of [5] for the formal dimension reduction of the plate model. A rigorous derivation for more general material models including the above St. Venant–Kirchhoff material can be found in [21] which leads to the same dimensionally reduced model up to a different prefactor in the energy density. For further details on admissible energies and relaxation effects of the dimension reduction, we refer the reader to [11, 10, 20, 21]. We remark that a more general, nonlinear dependence of the strain tensor on temperature can be used; for our experiments a linear relation turned out to be sufficient.

2.2. Bilayer configuration and evolution hypotheses

We consider the following configuration below. We let $\omega\subset\mathbb{R}^{2}$ be the flat parametric domain, $\delta>0$ be the thickness of the plate, and $\omega_{\delta}:=\omega\times(-\delta/2,\delta/2)$ be the plate in the undeformed configuration. We let $\mathbf{x}:=(\mathbf{x}^{\prime},x_{3})\in\omega_{\delta}$ denote a generic point and $t\in(0,\infty)$ time. We further indicate with $\mathbf{u}(\cdot,t):\omega_{\delta}\to\mathbb{R}^{3}$ the deformation of the plate and $\gamma_{\delta}(t)=\mathbf{u}(\cdot,t)(\omega_{\delta})$ the deformed configuration of the plate at time $t$.

For the evolution model considered below we assume that inertial effects are insignificant and that no mechanical dissipation occurs so that the plate immediately adjusts to temperature changes. This yields a dynamics driven by temperature, which mathematically entails that $\Theta$ solves a suitable diffusion equation on the flexible surface and $\mathbf{u}(\cdot,t)$ is a minimizer of an elastic energy associated with the density $W_{\Theta}(\nabla\mathbf{u}(\cdot,t))$. Due to the assumptions of thermal insulation of the bilayer top and bottom and of mid-surface inextensibility, it turns out that the diffusion equation decouples from the deformation equation. We examine these two models below and outline reduced models for the vanishing thickness limit $\delta\to 0$.

2.3. Reduced bilayer plate model

We adjust the simplified energy (1) for hyperelastic materials to model thin bilayers and derive effective energies describing large deformations. In particular, we consider two layers of materials glued on top of each other with different thermal material constants such as the thermal expansion coefficients. We assume that one material expands and the other compresses above a critical temperature, which we assume to be zero for simplicity. Therefore, we consider (1) with $F=\nabla\mathbf{u}$

[TABLE]

where the coefficients $\alpha(\mathbf{x}),\mu(\mathbf{x})$ jump across $x_{3}=0$

[TABLE]

We follow [11, 22, 5] to identify a dimensionally reduced model corresponding to the limit $\delta\to 0$. Formally, this is based on the Kirchhoff assumption that the actual deformation $\mathbf{u}$, subject to given forces and boundary conditions can be, up to higher order contributions, represented as

[TABLE]

with a mapping $\mathbf{y}:\omega\times(0,\infty)\to\mathbb{R}^{3}$ that describes the deformation of the midplane $\omega\times\{0\}$ and a vector field $\mathbf{b}:\omega\times(0,\infty)\to\mathbb{R}^{3}$ that is normal to the deformed midplane $\gamma(t)=\mathbf{y}(\omega,t)$. This means that fibers perpendicular to $\omega$ in the undeformed configuration remain normal to the mid-surface $\gamma(t)$ in the deformed configuration. For ease of presentation we do not write the argument $t$ explicitly in the remainder of this subsection. Inserting the corresponding deformation gradient

[TABLE]

into the scaled elastic energy functional

[TABLE]

leads to

[TABLE]

The scaling of the elastic energy by $\delta^{-3}$ corresponds to deformations that describe a bending behavior of the thin plate and is crucial for identifying limiting equations. We assume that $I_{\Theta}[\mathbf{u}]$ remains bounded as $\delta\to 0$, carry out the integration in $x_{3}$ direction, and deduce necessary scaling properties of terms arising in the energy functional; details can be found in [5]. The first necessary condition to have a finite limit as $\delta\to 0$ is that $\mathbf{b}$ has length $|\mathbf{b}|=1\pm\alpha\Theta$ in the upper and lower layers. Assuming that $\alpha/\delta$ is finite, which is explained below, shows that for $\delta\to 0$ we have that $\mathbf{b}$ equals the unit normal ${\boldsymbol{\nu}}$ to the deformed midplane $\gamma$ and that (\nabla^{\prime}\mathbf{b})^{\rotatebox[origin={c}]{180.0}{\footnotesize\perp}}\mathbf{b}=0. Let $G$ and $H$ be the first and second fundamental forms of $\gamma$, i.e., the symmetric matrices

[TABLE]

Then we find the second necessary condition $G=I_{2}$ for a finite limit $\delta\to 0$,

[TABLE]

i.e., $\mathbf{y}:\omega\to\mathbb{R}^{3}$ defines an isometric deformation of the midplane $\omega$. With this we derive the central identity

[TABLE]

in which $\overline{\mu}$ is the average of $\mu(\cdot,\pm x_{3})$, $\theta(\mathbf{x}^{\prime})=\Theta(\mathbf{x}^{\prime},0)$ is the effective temperature that obeys a reduced diffusion equation discussed below, and $\overline{\alpha}$ is the effective thermal expansion coefficient per unit thickness given by

[TABLE]

Note that $\overline{\alpha}\theta$ has the unit of a curvature, i.e. $\text{\,}{\mathrm{m}}^{-1}$.

We thus infer that the reduced energy functional takes the form

[TABLE]

up to constant terms that we ignore from now on. Assuming that $\alpha$ is comparable with the plate thickness is necessary to avoid delamination and obtain a total bending effect of the thin bilayer plate resulting from a compressive and expansive behavior of the individual layers. This also shows that small material differences can lead to large deformations.

2.4. Reduced diffusion model

As in the mechanical part we allow the two different materials in the upper and lower layers of the plate to possess different physical constants, namely heat capacities $\sigma^{\pm}$ and conductivities $\kappa^{\pm}$. We define

[TABLE]

to be their mean values in the transversal direction. We make the key simplifying assumptions that the plate is thermally insulated on top and bottom and that diffusion is taking place mostly in tangential direction of the plate. This means that the gradient of temperature, and so the heat flux $\mathbf{q}=-\kappa\nabla\Theta$ which obeys the Fourier law of diffusion, satisfies the orthogonal decomposition

[TABLE]

where for $(\mathbf{x}^{\prime},x_{3})\in\omega_{\delta}$

[TABLE]

as well as

[TABLE]

We start with the energy balance in the slender 3d set $\gamma_{\delta}(t):=\mathbf{u}(\cdot,t)(\omega_{\delta})$, i.e., in the deformed plate with positive thickness $\delta$ where $\mathbf{u}$ is defined in (2):

[TABLE]

Using Reynolds’ transport theorem, cf., e.g., [9], we can write the left-hand side as follows:

[TABLE]

Using (6), we obtain for the second summand in the right-hand side as $\delta\to 0$

[TABLE]

where $\theta(\mathbf{x}^{\prime},t)=\Theta(\mathbf{x}^{\prime},0,t)$ is the temperature in the midsurface $\gamma(t)$ and $\overline{\sigma}$ is defined in (5). It remains to evaluate the term $\operatorname{div}\partial_{t}\mathbf{u}$.

To do so, we set $\mathbf{F}(\mathbf{x})=\mathbf{f}(\mathbf{u})$ for a vector-valued function $\mathbf{f}:\gamma_{\delta}\to\mathbb{R}$ and use (2) to deduce

[TABLE]

Since the divergence is the trace of the gradient, we invoke the decomposition (6) for each component of $\mathbf{f}$ together with the preceding expression to obtain

[TABLE]

To apply this formula to $\mathbf{f}=\partial_{t}\mathbf{u}=\partial_{t}\mathbf{y}+x_{3}\,\partial_{t}\mathbf{b}$ we realize that $\partial_{x_{3}}\partial_{t}\mathbf{u}=\partial_{t}\mathbf{b}$ whence $\mathbf{b}\cdot\partial_{x_{3}}\partial_{t}\mathbf{u}=\mathbf{b}\cdot\partial_{t}\mathbf{b}=\frac{1}{2}\partial_{t}|\mathbf{b}|^{2}=0$ and

[TABLE]

Therefore, as $\delta\to 0$ we arrive at

[TABLE]

We now deal with the flux term $\mathbf{q}=-\kappa\nabla\Theta$. We first write

[TABLE]

where $\operatorname{div}_{\gamma}^{\perp}\mathbf{q}={\boldsymbol{\nu}}\cdot\nabla\mathbf{q}\,{\boldsymbol{\nu}}$ in view of (9). Integrating by parts in the normal direction ${\boldsymbol{\nu}}$, and using the vanishing Neumman boundary conditions $\mathbf{q}^{+}\cdot{\boldsymbol{\nu}}=\mathbf{q}^{-}\cdot{\boldsymbol{\nu}}=0$ on top and bottom of the plate $\gamma_{\delta}$ gives

[TABLE]

On the other hand, in light of (6) and (7) we deduce as $\delta\to 0$

[TABLE]

where $\overline{\kappa}$ is defined in (5) because

[TABLE]

Collecting all the previous results, we conclude that the limit of (8) as $\delta\to 0$ is the surface conservation equation

[TABLE]

This equation is consistent with the diffusion equation on a surface derived in [9].

We conclude with a simple extension of (11) which accounts for diffusion transversal to the plate. Suppose that the normal fluxes $\mathbf{q}^{+},\mathbf{q}^{-}$ on top and bottom of the plate $\gamma_{\delta}(t)$ do not vanish but rather scale proportional to $\delta$. In this case, (10) reduces to

[TABLE]

In other works, the function $f$ acts as an effective source term in the surface diffusion equation. In addition, the following limit does no longer vanish

[TABLE]

where $V={\boldsymbol{\nu}}\cdot\partial_{t}\mathbf{y}$ is the normal velocity and $h=-{\rm tr}H$ is the mean curvature of $\gamma(t)$. This leads to the following variant of (11)

[TABLE]

which of course must be supplemented with a diffusion equation in the surroundings of $\gamma(t)$ to determine the quantity $\partial_{\boldsymbol{\nu}}\theta$. In contrast to (11), this form of the PDE does couple diffusion and plate geometry.

3. Governing PDEs: Weak forms

The governing equations are formulated in the cylinder $\omega\times(0,T)$ and the independent variables are denoted by $(\mathbf{x},t)\in\omega\times(0,T)$ for simplicity.

3.1. Plate equation

We first derive the Euler-Lagrange equation for a minimizer of the bending energy (4) given a fixed temperature distribution.

In (4), we omit writing the variable $\mathbf{x}$ and write

[TABLE]

For isometries, the $ij$-th element of $H$, namely $h_{ij}=\partial_{i}\partial_{j}\mathbf{y}\cdot{\boldsymbol{\nu}}$, satisfies the key relation

[TABLE]

hence $\partial_{i}\partial_{j}\mathbf{y}$ is parallel to ${\boldsymbol{\nu}}$. This immediately implies equality of the Frobenius norm of the second fundamental form $H$ and the Hessian of $\mathbf{y}$

[TABLE]

Developing the square yields the following equivalent expression for $E[\mathbf{y}]$

[TABLE]

Since ${\boldsymbol{\nu}}=\partial_{1}\mathbf{y}\times\partial_{2}\mathbf{y}$ for isometries, in view of (12) we deduce the expression

[TABLE]

Consequently, $\mathbf{y}$ is a minimizer of the energy

[TABLE]

subject to the isometry constraint

[TABLE]

Computing the first variation of these two relations with respect to $\mathbf{y}$ yields the weak form of the Euler-Lagrange equation

[TABLE]

for any test function $\mathbf{w}$ satisfying the linearized isometry condition

[TABLE]

The identity (14) can be simplified further because the third term on the left-hand side vanishes. To see this, we first recall from (12) that second derivatives $\partial_{i}\partial_{j}\mathbf{y}$ of $\mathbf{y}$ are parallel to ${\boldsymbol{\nu}}=\partial_{1}\mathbf{y}\times\partial_{2}\mathbf{y}$ for isometries, whence the mean curvature $h=-{\rm tr}(HG^{-1})=-{\rm tr}H$ satisfies

[TABLE]

We next observe that $\partial_{i}\mathbf{y}\cdot\partial_{j}\mathbf{y}=0$ for $i\neq j$ according to (13) and $\partial_{i}\mathbf{y}\cdot\partial_{i}\mathbf{w}=0$ because of (15). Using now the formula $(a\times b)\cdot(c\times d)=(a\cdot c)(b\cdot d)-(a\cdot d)(b\cdot d)$ yields

[TABLE]

This implies that $\mathbf{y}$ is a solution of the simplified fourth order equation

[TABLE]

for all test functions $\mathbf{w}$ satisfying the linearized isometry condition (15).

3.2. Diffusion equation

We now turn to the diffusion equation (11). Since its derivation is valid for any patch in $\gamma(t)$, we realize that it implies the strong form

[TABLE]

where $D_{t}\theta=\partial_{t}\theta+\nabla_{\gamma}\theta\cdot\partial_{t}\mathbf{y}$ is the material derivative of $\theta$, cf. [9]. We split the boundary $\partial\gamma(t)$ of $\gamma(t)$ in two disjoint pieces $\partial_{D}\gamma(t)$ and $\partial_{R}\gamma(t)$, where $\partial_{D}\gamma(t)$ is the portion of $\partial\gamma(t)$ where we prescribe the temperature and $\partial_{R}\gamma(t)$ is where we impose the Robin condition

[TABLE]

Multiplying (16) by any test function $\varphi$ with vanishing material derivative $D_{t}\varphi=0$ and vanishing on $\partial_{D}\gamma(t)$, and integrating over $\gamma(t)$ leads to

[TABLE]

We thus have that

[TABLE]

this form of the diffusion equation on $\gamma(t)$ is due to [9]. In the present context, this equation simplifies further because $\mathbf{y}$ is an isometry for all $t$. In fact, since the first fundamental form satisfies G=[\nabla\mathbf{y}]^{\rotatebox[origin={c}]{180.0}{\footnotesize\perp}}\nabla\mathbf{y}=I_{2} we can write with $\mathbf{z}=\mathbf{y}(\mathbf{x})$

[TABLE]

and

[TABLE]

because $\nabla_{\gamma}\varphi(\mathbf{z},\cdot)=\nabla\mathbf{y}G^{-1}\nabla\varphi(\mathbf{y}(\mathbf{x},\cdot),\cdot)$. This allows us to express (17) in the parametric domain $\omega$ and avoid the dependence on $t$ in the test function $\varphi$. We thus get the following simple weak form of the diffusion equation

[TABLE]

i.e., it is sufficient to solve the diffusion equation in the reference configuration $\omega$. Moreover, the diffusion equation decouples from the plate equation due to the isometry property and the vanishing Neumann boundary conditions for temperature assumed on top and bottom of the plate.

4. Numerical scheme

4.1. Time discretization

We next describe a discrete-time scheme to compute the evolution of thermally induced bilayer bending effects. We consider clamped boundary conditions for the mechanical equation and Dirichlet conditions for the diffusion process on the subset $\partial_{D}\omega\subset\partial\omega$. A Robin-type boundary condition is imposed on the remaining part $\partial_{R}\omega=\partial\omega\setminus\partial_{D}\omega$ while no explicit boundary conditions are imposed on the deformation on this part. We further impose the condition that deformations are contained in a convex set $K$ which models the presence of an obstacle, e.g., $K=\{\mathbf{y}=(y_{i})_{i=1}^{3}\in L^{2}(\omega;\mathbb{R}^{3}):y_{3}\leq 1\}$, or $K$ is the entire space $K=L^{2}(\omega;\mathbb{R}^{3})$ in the case of no obstacle. The $L^{2}$ scalar product on $\omega$ of functions or vector fields is denoted by $(\cdot,\cdot)$. If the $L^{2}$ scalar product is taken on a set $A$ we write $(\cdot,\cdot)_{A}$. In addition, we let $\{\mathbf{e}_{1},\mathbf{e}_{2}\}$ be the canonical unit vectors in $\omega$ and use the notation $d_{t}\theta^{k+1}$ to indicate the scaled backward difference

[TABLE]

Algorithm 1** (abstract time stepping).**

*Let $\tau>0$ be a uniform time step, set $\theta^{0}:=0$, $\mathbf{y}^{0}:={\rm id}$, $\nabla\mathbf{y}^{0}:=[\mathbf{e}_{1},\mathbf{e}_{2}]$ and $k:=0$.

(1) Compute $\theta^{k+1}\in H^{1}(\omega)$ with $\theta^{k+1}|_{\partial_{D}\omega}=\theta_{D}$ and*

[TABLE]

*for all $\varphi\in H^{1}_{D}(\omega):=\{\psi\in H^{1}(\omega):\psi|_{\partial_{D}\omega}=0\}$.

(2) Compute a minimizer $\mathbf{y}^{k+1}\in H^{2}(\omega;\mathbb{R}^{3})$ of the functional*

[TABLE]

subject to the obstacle constraint $\mathbf{y}\in K$, the boundary conditions $\mathbf{y}^{k+1}=[{\rm id},0]$ and $\nabla\mathbf{y}^{k+1}=\nabla[\mathbf{e}_{1},\mathbf{e}_{2}]$ on $\partial_{D}\omega$ and the linearized isometry condition

[TABLE]

(3) Increase $k\to k+1$ and continue with (1).

Note that the operators $\nabla$ and $D^{2}$ entail partial derivatives with respect to the parametric variables $\mathbf{x}=(x_{1},x_{2})\in\omega$ only. To obtain a practical, semi-implicit version of Algorithm 1 we use the identities derived in Section 3. Given an approximation $\mathbf{y}^{k}$ we define a tangent space relative to the isometry constraint and the boundary conditions by

[TABLE]

i.e., vector fields $\mathbf{v}$ satisfying homogeneous clamped boundary conditions and the linearized isometry constraint (15). Note that the space $H^{2}_{D}(\omega)$ consists of all functions that vanish along with their gradients on $\partial_{D}\omega$.

Obstacle constraint. We deal with the constraint $\mathbf{y}\in K$ via a variable splitting and penalization of such splitting in the $L^{2}$-norm. We thus introduce the auxiliary variable $\mathbf{s}\approx\mathbf{y}$, penalize the deviation of $\mathbf{s}$ from $\mathbf{y}$ by adding the penalty $L^{2}$-term $\frac{1}{12\varepsilon}\|\mathbf{y}-\mathbf{s}\|^{2}$ to the energy $I[\mathbf{y}]$, i.e. consider

[TABLE]

and minimize $J[\mathbf{y},\mathbf{s}]$ separately over $\mathbf{y}$ and $\mathbf{s}$ imposing that $\mathbf{s}\in K$. Note that minimizing with respect to $\mathbf{s}$ leads to the $L^{2}(\omega)$-orthogonal projection of $\mathbf{y}$ onto $K$, which we denote by $\mathbf{s}=\Pi_{K}(\mathbf{y})$, whereas $\mathbf{y}$ is unconstrained. We assume that the undeformed plate does not intersect the obstacle.

Algorithm 2** (practical time stepping).**

*Let $\tau>0$, set $\theta^{0}:=0$, $\mathbf{y}^{0}:=[{\rm id},0]$, $\nabla\mathbf{y}^{0}:=[\mathbf{e}_{1},\mathbf{e}_{2}]$, $\mathbf{s}^{0}:=\mathbf{y}^{0}$, and $k:=0$.

(1) Compute $\theta^{k+1}\in H^{1}(\omega)$ with $\theta^{k+1}|_{\partial_{D}\omega}=\theta_{D}$ and*

[TABLE]

*for all $\varphi\in H^{1}_{D}(\omega)$.

(2) Compute $\mathbf{v}^{k+1}\in\mathcal{F}[\mathbf{y}^{k}]$ such that*

[TABLE]

for all $\mathbf{w}\in\mathcal{F}[\mathbf{y}^{k}]$ and set $\mathbf{y}^{k+1}:=\mathbf{y}^{k}+\tau\mathbf{v}^{k+1}$. Set

[TABLE]

(3) Increase $k\to k+1$ and continue with (1).

The precise stability properties of Algorithm 2 appear difficult to identify. In particular, due to the lack of control on a discrete time-derivative for $(\mathbf{y}^{k})$, the inconsistencies related to the explicit treatment of some terms cannot be controlled directly. If we write $\mathbf{y}=\mathbf{y}^{k}+\tau\mathbf{v}$ with $\mathbf{v}\in\mathcal{F}[\mathbf{y}^{k}]$, then the equation in Step (2) is the Euler-Lagrange equation of the energy

[TABLE]

without obstacle constraint on $\mathbf{y}$. In contrast, $\mathbf{s}^{k+1}$ is the $L^{2}$-projection of $\mathbf{y}^{k+1}$ onto $K$, namely $\mathbf{s}^{k+1}=\Pi_{K}(\mathbf{y}^{k+1})$, which does not involve any solve. The decoupling of $\mathbf{y}^{k+1}$ and $\mathbf{s}^{k+1}$ in Step (2) is motivated by separate convexity properties of $J[\mathbf{y},\mathbf{s}]$ in each argument. A simultaneous minimization of $J[\mathbf{y},\mathbf{s}]$ in $\mathbf{y}$ and $\mathbf{s}$ can be iteratively realized by repeating the two substeps in Step (2) where the term $\varepsilon^{-1}(\mathbf{s}^{k},\mathbf{w})$ is repeatedly replaced by $\varepsilon^{-1}(\mathbf{s}^{k+1},\mathbf{w})$ until the subiteration becomes stationary.

4.2. Space discretization

Let $\mathcal{T}_{h}$ be a partition of the reference domain $\omega$ into quadrilaterals with diameters comparable to $h$. We use the standard lowest order $H^{1}$ conforming finite element space $\mathbb{Q}_{1}$ to discretize the diffusion equation. We solve approximately the fourth order nonlinear bending problems with $H^{2}$ nonconforming discrete Kirchhoff quadrilaterals [5]. The key idea in their construction is the use of two $H^{1}$ conforming finite element spaces to approximate deformations and deformations gradients together with a discrete (or reduced) gradient operator $\nabla_{h}$ that connects the spaces. We adopt the description of the element from [5] which is motivated by the triangular version considered in [3, 7, 8]. We let $\mathbb{Q}_{r}(T)$ and $\mathbb{P}_{r}(T)$ denote the set of polynomials on $T\in\mathcal{T}_{h}$ of partial degree $r$ on each variable and of total degree $r$, respectively. Let $\mathcal{N}_{h}$ be the set of vertices of elements in $\mathcal{T}_{h}$ and $\mathcal{E}_{h}$ be the set of edges in $\mathcal{T}_{h}$. For every $E\in\mathcal{E}_{h}$ we let $\mathbf{n}_{E}$ be a unit normal vector to $E$ and $\mathbf{z}_{E}$ be the midpoint of $E$.

Definition 1** (discrete spaces and operators).**

(i) Define the discrete spaces

[TABLE]

(ii) Let $\widehat{\mathcal{I}}_{h}^{2}:[H^{2}(\omega)]^{2}\to\mathbb{G}_{h}$ be the interpolation operator defined by

[TABLE]

*The operator $\widehat{\mathcal{I}}_{h}^{2}$ is also well-defined for discrete vector fields $\psi\in\nabla\mathbb{W}_{h}$.

(iii) Let $\nabla_{h}:H^{3}(\omega)\to\mathbb{G}_{h}$ be the discrete gradient operator defined by*

[TABLE]

The operator $\nabla_{h}$ is also well-defined for discrete functions $w\in\mathbb{W}_{h}$.

Note that $\nabla\mathbb{W}_{h}$ is a space of discontinuous vector fields containing $\mathbb{Q}_{2}(T)$ for all $T\in\mathcal{T}_{h}$, whereas $\nabla_{h}\mathbb{W}_{h}$ is a smaller conforming space.

For an efficient numerical treatment of nonlinearities such as the projection operator $\Pi_{K}$ we define the discrete inner product $(\cdot,\cdot)_{h}$ for piecewise continuous functions $\phi,\psi\in\Pi_{T\in\mathcal{T}_{h}}[C^{0}(T)]^{\ell}$

[TABLE]

To define our numerical scheme we also define

[TABLE]

The discrete counterpart of the isometry condition (13) is imposed at the nodes $\mathcal{N}_{h}$ of the mesh $\mathcal{T}_{h}$, which leads to the following discrete linearization:

[TABLE]

The following scheme then only requires the solution of linear systems of equations and the evaluation of low dimensional projections.

Algorithm 3** (fully practical scheme).**

*Let $\tau>0$, set $\theta_{h}^{0}:=0$, $\mathbf{y}_{h}^{0}:=[{\rm id},0]$, $\nabla\mathbf{y}_{h}^{0}:=[\mathbf{e}_{1},\mathbf{e}_{2}]$, $\mathbf{s}_{h}^{0}:=\mathbf{y}_{h}^{0}$, and $k:=0$.

(1) Compute $\theta_{h}^{k+1}\in\mathbb{V}_{h}$ with $\theta_{h}^{k+1}|_{\partial_{D}\omega}=\theta_{D}$ and*

[TABLE]

*for all $\varphi_{h}\in\mathbb{V}_{h}\cap H^{1}_{D}(\omega)$.

(2) Compute $\mathbf{v}_{h}^{k+1}\in\mathcal{F}_{h}[\mathbf{y}_{h}^{k}]$ such that*

[TABLE]

for all $\mathbf{w}_{h}\in\mathcal{F}_{h}[\mathbf{y}_{h}^{k}]$ and set $\mathbf{y}_{h}^{k+1}:=\mathbf{y}_{h}^{k}+\tau\mathbf{v}_{h}^{k+1}$. Set

[TABLE]

*for all $\mathbf{z}\in\mathcal{N}_{h}$.

(3) Increase $k\to k+1$ and continue with (1).*

The definitions of both $\mathcal{F}_{h}[\mathbf{y}_{h}^{k}]$ and $\mathbf{y}_{h}^{k+1}$ imply that $\mathbf{y}_{h}^{k+1}=[{\rm id},0],\nabla y_{h}^{k+1}=[\mathbf{e}_{1},\mathbf{e}_{2}]$ on $\partial_{D}\omega$ for all $k\geq 0$.

5. Numerical experiments

We aim at numerically investigating practically relevant scenarios for our model problem. For this purpose we work with physical units and realistic material parameters in the following. These apply to polydimethylsiloxane, polyvinyl-alcohol, polystyrene, soft polyvinylchloride, and soft polyurethane materials.

5.1. Material parameters

The material parameters involved in our mathematical model belong to the following ranges for typical polymer materials:

[TABLE]

The heat capacity is then $\sigma=\varrho c_{v}$. Except for the thermal expansion coefficient, which has opposite signs in the upper and lower layer, we use the same material parameters for the two layers. Unless stated otherwise we use

[TABLE]

We consider different geometries $\omega$ whose diameters are on the order of a few milimeters but with the same thickness

[TABLE]

The resulting effective material parameters for the reduced model are then

[TABLE]

Here, we used (3), (5), $\overline{\eta}=(\eta^{+}+\eta^{-})/2$, and the formula

[TABLE]

from [10] which results from a more general dimension reduction.

5.2. Discretization parameters

The choice of meshsize $h$ is dictated by space resolution, which improves with the use of mesh refinement in regions of major activity, typically near the hinges. The choice of time step $\tau$ and penalty parameters $\varepsilon$ is trickier to yield to realistic physical rather than numerical effects. To see why, we present now a nondimensional analysis. We choose characteristic length $\ell$, temperature $\theta_{0}$, bending coefficient $\mu_{0}$ and time $T$, and denote the scaled variables

[TABLE]

and scaled parameters

[TABLE]

We next rewrite the funcional in (18) in terms of the new variables:

[TABLE]

This reveals that the effective penalty parameter is

[TABLE]

Moreover, since the evolution is dictated by the diffusion equation, we would like the elastic energy minimization to reflect the quasi-stationary nature of the plate evolution. If the $k$-iterate $\widehat{\mathbf{y}}^{k}$ is far from the obstacle, then $\widehat{\mathbf{s}}^{k}=\Pi_{K}(\widehat{\mathbf{y}}^{k})=\widehat{\mathbf{y}}^{k}$ and we can rewrite the penalty term as follows

[TABLE]

where d_{\widehat{t}}\,\widehat{\mathbf{y}}^{k+1}=\widehat{\tau}^{-1}\big{(}\widehat{\mathbf{y}}^{k+1}-\widehat{\mathbf{y}}^{k}\big{)} is the discrete time derivative. Consequently, to capture the expected physical behavior we impose the condition

[TABLE]

In other words, the two limits $\tau\to 0$ and $\varepsilon\to 0$ do not commute and the former has to take place before the latter for model consistency.

5.3. Bilayer micro-scale valves and switching devices

Adaptive bilayer materials are promising for micro-valves for fluid flow management systems, where flapping or flexure mechanisms are used to control the flow. For example, thermally actuated micro-valves made of bilayers [12] have been proposed where closing and opening a valve aperture is controlled through heating. Another example is using piezoelectric ceramics for micro-valves [16]. As ceramics is brittle, limited deformation can be obtained in closing and opening the valves. Alternatively, polymer bilayers can be actuated to attain relatively large deformations and used in micro-valves.

As proof of concept, we propose in this section a numerical study of bilayer switching devices consisting of bilayer hinges assembled with flexible plates. Such thermally operated devices have successfully been constructed as well. Upon activation the bilayer rotates the plate in order to touch an object. We model this by considering the square $\omega=(-1,1)^{2}$ with side lengths $2.0\text{\,}\mathrm{mm}$ which is composed of a flexible single layer and a thin bilayer strip of width $\pi/40$, cf. Figure 2. The edge $\partial_{D}\omega=\{-1\}\times[-1,1]$ of the bilayer strip which is not connected to the plate is assumed to be clamped, i.e., the deformation $\mathbf{y}$ is constrained to satisfy the boundary condition

[TABLE]

The plate is initially at the critical temperature, i.e., we impose the condition

[TABLE]

The bending mechanism is triggered by heating the plate from the clamped side via the prescribed temperature

[TABLE]

The hinge thickness is so small that heat diffuses quite fast. To prevent a very rapid motion of the plate, which reacts instantaneously to the hinge bending, the factor $\min\left(1,\frac{t}{5}\right)$ ensures a smooth and slow transition from the initial temperature $\theta=$0\text{,}\mathrm{\SIUnitSymbolCelsius} to the desired one of $100\text{\,}\mathrm{\SIUnitSymbolCelsius}$. We assume that the plate is thermally insulated and that the deformation is free on the remaining part $\partial_{R}\omega:=\partial\omega\setminus\partial_{D}\omega$ of the boundary, i.e., we use a vanishing effective heat transfer coefficient $\overline{\eta}$. The plate rigidity is uniform throughout the plate and characterized by a bending coefficient $\overline{\mu}=$2.0\text{\times}{10}^{3}\text{\,}{\rm MPa}. A flat obstacle that models the contact object is modeled by the set

[TABLE]

and the projection $\Pi_{K}$ onto $K$ is computed nodewise for $\mathbf{y}_{h}=(y_{i})_{i=1}^{3}$:

[TABLE]

For our numerical experiments we use a partition of $\omega$ that results from 6 uniform quad refinements of $\omega$. We use uniform time steps $\tau=$3.0\text{\times}{10}^{-3}\text{,}\mathrm{s} and penalization parameter $\varepsilon=$4\text{\times}{10}^{-6}\text{\,}\mathrm{m}\mathrm{m}^{4}{\rm MPa}^{-1}. Since $\mu_{0}=2\times 10^{3}{\rm MPa}$ and $T=10s$, this choice is consistent with (19), namely $\tau^{2}\ll 0.5\times 10^{-1}\,s^{2}$. Figure 2 shows the evolution from the initial configuration.

To assess the intricate influence of the penalty parameter $\varepsilon$, we consider the same setting except for values of $\varepsilon=4\times 10^{-j}\,$\text{,}\mathrm{m}\mathrm{m}^{4}{\rm MPa}^{-1}$$ for $j=4,...,9$. Cuts along the plane $\{x_{2}=0\}$ (perpendicular to the obstacle and the clamped side) of the plates stationary deformations are displayed in Figure 3. Plates are considered in a stationary state whenever

[TABLE]

On the one hand, the smallest value $\varepsilon$ corresponding to $j=9$ results in a plate which does not cross the obstacle but is influenced by the presence of the obstacle already when still far from the plate. On the other hand the largest value of $\varepsilon$ corresponding to $j=6$ results in a effective deformations without obstacle. Values of $\varepsilon$ for parameters $j=7,8$ yield a deformation crossing the planar obstacle by no more that the finite element meshsize.

It is worth mentioning that except the case $j=9$, all the stationary states (defined according to (21)) are reached at $T=$136.67\text{,}\mathrm{s} (44689 time iterations) while $j=9$ required $T=$182.739\text{\,}\mathrm{s} (60913 time iterations). However, the dynamics for smaller $\varepsilon$ is influenced by the obstacle at early stages as illustrated in Figure 4. This is consistent with (19).

5.4. Dog-ear formation

To explore possible failure of controlled production of microtubes reported in [28, 33, 34] we consider a bilayer square $\omega=(-1,1)^{2}$ of side lengths $2.0\text{\,}\mathrm{mm}$ that is clamped and thermally insulated on the side $\partial_{D}\omega=\{-1\}\times[-1,1]$, namely $\mathbf{y}_{h}$ satisfies (20) and $\partial_{\mathbf{n}}\theta=0$ on $\partial_{D}\omega$. On the remaining part $\partial_{R}\omega=\partial\omega\setminus\partial_{D}\omega$ the plate is free and heated via external heat transfer obeying Newton’s law of cooling, i.e., the normal flux is proportional to temperature difference

[TABLE]

relative to the ambient temperature $\theta_{ext}=$50.0\text{,}\mathrm{\SIUnitSymbolCelsius}$$. The experiment setup is sketched in Figure 5. For our numerical experiments we consider two different values of the effective thermal conductivity $\overline{\kappa}$ and heat capacity $\overline{\sigma}$ so that the resulting diffusivity differs by a factor 10, i.e., we consider

[TABLE]

We run Algorithm 3 (fully discrete scheme) for a uniform partitition of $\omega$ into subsquares with side lengths $h=2^{-6}\times$2.0\text{,}\mathrm{mm} and uniform time step $\tau=$5.0\text{\times}{10}^{-3}\text{\,}\mathrm{s}. Figure 5 displays snapshots of the approximate evolution for settings (a) and (b) for the respective times

[TABLE]

We observe that the smaller diffusivity case (a) leads to accumulation of heat and higher temperatures in the vicinity the boundary $\partial_{R}\omega$ and in particular at the two corners belonging only to $\partial_{R}\omega$. This results in a stronger localized bending behavior for (a) and the formation of dog-ears which then prevent a controlled rolling up into a cylindrical shape. For larger diffusivity (b), heat distributes faster within the plate and the bending behavior is less localized.

5.5. Self-assembling box

We consider an arrangement of rigid plates and bilayer hinges that realizes a self-assembling and self-opening box actuated by temperature. We refer the reader to [29, 13] for related experiments. The relevant practical applications on this type of folding systems are on deployable structures such as drug delivery systems, deployable shelters with prefabricated walls. As in Subsection 5.3, the plate has $\overline{\alpha}=0$ and $\overline{\mu}$ that is $20$ larger than that of the hinge, so as to mimic a rigid plate insensitive to temperature. The geometry is sketched in Figure 6.

The square rigid plates have side lengths $1.0\text{\,}\mathrm{mm}$ and the hinges have widths $\pi\text{/}48\text{\,}\mathrm{mm}$. The central square, which is connected to four other squares, is assumed to be permanently attached to a substrate, and a constant heat source of $75.0\text{\,}\mathrm{\SIUnitSymbolCelsius}\text{\,}{\mathrm{s}}^{-1}$ is applied on a circle of radius $0.25\text{\,}\mathrm{mm}$ centered on this square during a period of $19.0\text{\,}\mathrm{s}$. Note that in that case, the heat equation reads

[TABLE]

When temperature increases the bilayer hinges start bending while the rigid plates remain undeformed, cf. Figure 6. For sufficient heating a box is formed which can be unfolded via external cooling. For our simulation we chose a uniform time step $\tau=$0.5\text{,}\mathrm{s}$$ and a partition of $\omega$ into $3$ uniform quad refinements of each square plate, followed by $4$ quad refinements of elements in the hinges and additional local refinements to have at most one hanging node per side.

5.6. Deployable airfoil

We consider another arrangement of bilayer hinges and rigid plates shown in Figure 7 that models an initially flat device that can fold into an airfoil via external heating. Radical shape changes/folding of wing structures during flight can lead to high maneuverability of the plane and an efficient cruising that is adjustable to various flight envelopes. One example of morphing wing technology, i.e. z-shape folding, has been proposed by the Defense Advanced Research Projects Agency (DARPA), [17].

The physical parameters are the same as in Subsection 5.5 except for $\overline{\kappa}/\overline{\sigma}=$0.1\text{,}\mathrm{mm}\text{,}{\mathrm{}}^{2}\text{,}{\mathrm{s}}^{-1} and the middle plate is fixed. As indicated in the sketch of Figure [7](#S5.F7), the bilayer hinges bend into different directions, which is modeled by different signs of the effective expansion coefficients, i.e. $\overline{\alpha}=\pm 0.3$, or equivalently upon inverting the bilayer. The heat diffusion process is initiated by prescribing the temperature $\theta_{D}=$60.0\text{\,}\mathrm{\SIUnitSymbolCelsius} at the Dirichlet boundary parts $\partial_{D}\omega$ where the hinges meet $\partial\omega$, and vanishing Neumann condition on the rest $\partial_{R}\omega=\partial\omega\setminus\partial_{D}\omega$ (i.e. $\overline{\eta}=0$). A simulation of the thermally driven folding process is illustrated in Figure 7.

5.7. Particle encapsulation

Thermally controlled bilayers for transporting particles at microscales have been tested experimentally in [27]. Once deployed in a desirable place, enclosed particles may be released via external cooling of the device. Similar mechanisms can be triggered by changing the pH concentration of a surrounding liquid. This may find exciting and important medical applications in targeted drug delivery.

To simulate the essential effects we consider a star-shaped configuration of a bilayer as shown in Figure 8. To start the encapsulation process we use a Robin boundary condition on the entire boundary with external temperature

[TABLE]

The boundary is mechanically free but we fix the deformation at the vertices of one element in the mesh $\mathcal{T}_{h}$ which contains the midpoint of the center square to guarantee well-posedness of our method, i.e., uniquely defined iterates in Algorithm 3 (fully practical scheme). If the device center has the coordinate $(0,0)$, then the obstacle consists of $5$ spheres $B_{i}$, $i=0,...,4$, each of radius $0.24\text{\,}\mathrm{m}\mathrm{m}$ and centered at $(0.28,0.28,0.25)$, $(0.72,0.28,0.25)$, $(0.28,0.72,0.25)$, $(0.72,0.72,0.25)$ and $(0.5,0.5,0.5)$ (all in $$\text{,}\mathrm{m}\mathrm{m}$)$. At each vertex, the projection to the obstacle required by Algorithm 3 is approximated by the projection to the closest sphere. The obstacle models particles that are encapsulated by the deformed bilayer plate.

We construct the mesh $\mathcal{T}_{h}$ via 6 uniform quad refinements of the initial coarse partition of the domain $\omega$, indicated by the dashed lines in Figure 8. For our simulations we use the uniform time step $\tau=$2.5\text{\times}{10}^{-4}\text{,}\mathrm{s} and the penalization parameter $\varepsilon=$5.0\text{\times}{10}^{-8}\text{\,}\mathrm{m}\mathrm{m}^{4}{\rm MPa}^{-1}. Since $T=1s$ and $\mu_{0}=2\times 10^{3}{\rm MPa}$, this choice satisfies (19). Figure 8 depicts the encapsulation process of 5 non-penetrable and rigid spherical particles with radii $0.24\text{\,}\mathrm{mm}$. Penalization of the discrepancy between $y_{h}^{k+1}$ and $s_{h}^{k}$ in Algorithm 3 (fully practical scheme) does not prevent penetration of the obstacle, an effect that depends on the size of the penalization parameter $\varepsilon$ and the finite element meshsize and is discussed in Section 5.3. We refer to the lower part of Figure 8 that illustrates this feature. Reducing the size of $\varepsilon$ ameliorates this numerical artifact without making the algebraic system for $y_{h}^{k+1}$ of Algorithm 3 stiffer.

6. Conclusions

This paper develops a reduced model for the thermal actuation of bilayer plates consisting of a 4-th order nonlinear PDE and a 2-nd order diffusion in the mid-surface of a thin plate. The evolution is dictated by the latter whereas the former corresponds to a quasi-stationary reaction of the plate. We design a novel and effective finite element method for simulating large 3D deformations of slender planar structures, including the presence of obstacles. Several simulations illustrate the virtues of this approach.

The planar structures must be macroscopically compliant to allow for 3D shape changes with relatively small external stimuli. We study thin polymer materials composed of two layers which respond differently to thermal actuation, resulting in out of plane curvature changes and extremely large 3D deformations. The thin layer structures experience relatively small strain/stretch while undergoing large rotation. The interface between the two layers, or mid surface, is thus assumed to be perfectly bonded and inextensible, which leads to isometric deformations. The equations governing heat conduction and mechanical deformation decouple, i.e., the thermal effect induces plate deformations, but the latter do not affect heat conduction. This simplifies the solution process.

The algorithm takes advantage of several geometric properties valid for isometries, and advances in time with a semi-implicit scheme; the algebraic equations to be solved at each time step are linear. The only restrictions among discretization parameters come from achieving physically realistic dynamics in the presence of obstacles. Simulations of significant practical applications examine the feasibility of certain geometries and material parameters of polymers in attaining dramatic shape changes. This methodology could be used at the design stage, upon choosing material and geometric parameters as well as external stimuli and simulating shape reconfigurations prior to fabricating structures. This methodology is general and can be extented for more complex material models for active materials.