Numerical approximation of von K\'{a}rm\'{a}n viscoelastic plates

TL;DR

This paper develops a theoretical framework for approximating viscoelastic von Kármán plates using finite element methods, proving convergence of discretizations and demonstrating results through computational experiments.

Contribution

It introduces an abstract convergence result for discretizations of metric gradient flows and applies it to finite element approximation of viscoelastic plates.

Findings

Convergence of discretized schemes to curves of maximal slope.

Successful finite element approximation of viscoelastic von Kármán plates.

Computational experiments validate theoretical results.

Abstract

We consider metric gradient flows and their discretizations in time and space. We prove an abstract convergence result for time-space discretizations and identify their limits as curves of maximal slope. As an application, we consider a finite element approximation of a quasistatic evolution for viscoelastic von K\'{a}rm\'{a}n plates. Computational experiments are provided, too.

Equations

$$\displaystyle-{\rm div}\Big{(}\partial_{F}W(\nabla y)+\partial_{\dot{F}}R(\nabla y,\partial_{t}\nabla y)\Big{)}=f\text{ in $[0,T]\times\Omega$.}$$

$$R(F,\dot{F})=\tilde{R}(C,\dot{C})$$

$$\displaystyle\nabla_{h}y:=\Big{(}\nabla^{\prime}y,\frac{1}{h}y_{,3}\Big{)}=\nabla w.$$

$$({\nabla }_h^{2}y{)}_{ijk}:=h^{-{\delta }_{3j}-{\delta }_{3k}}({\nabla }^{2}y{)}_{ijk}=({\nabla }^{2}w{)}_{ijk}={\partial }_{jk}^2{w}_i\text{for i,j,k∈{1,2,3}},$$

$$\displaystyle\begin{split}(i)&\ \ W\text{ continuous and $C^{3}$ \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}in a neighborhood of $SO(3)$},\\ (ii)&\ \ \text{frame indifference: }W(QF)=W(F)\text{ for all }F\in\mathbb{R}^{3\times 3},Q\in SO(3),\\ (iii)&\ \ W(F)\geq c\operatorname{dist}^{2}(F,SO(3)),\ W(F)=0\text{ iff }F\in SO(3),\end{split}$$

$$\begin{array}{rl}(i)& \text{frame indifference: }P(QZ)=P(Z)\text{ for all }Z\in {\mathbb{R}}^{3\times 3\times 3},Q\in SO(3),\\ (ii)& \text{P is convex and C1},\\ (iii)& \text{growth condition: For all Z∈R3×3×3 we have }\\ & c_1|Z{|}^p\le P(Z)\le c_2|Z{|}^p,|{\partial }_{Z}P(Z)|\le c_2|Z{|}^{p-1}\end{array}$$

$$R(F,\dot{F})=R(QF,Q(\dot{F}+AF))\forall Q\in SO(3),A\in {\mathbb{R}}_{\mathrm{skew}}^{3\times 3}$$

$$R(F,\dot{F}):=\underset{\epsilon \to 0}{\lim }\frac{1}{2{\epsilon }^2}{D}^2(F+\epsilon \dot{F},F)=\frac{1}{4}{\partial }_{F_1^2}^2{D}^2(F,F)[\dot{F},\dot{F}]$$

$$(i)$$

$$(ii)$$

$$(iii)$$

$$(iv)$$

$$(v)$$

$$\forall Q_1,Q_2\in SO(3),\forall F_1,F_2\in G{L}_{+}(3),$$

$$(vi)$$

$$\displaystyle\mathfrak{S}_{h}=\Big{\{}y\in W^{2,p}(\Omega;\mathbb{R}^{3}):\ y(x^{\prime},x_{3})=\begin{pmatrix}x^{\prime}\\ hx_{3}\end{pmatrix}\text{for }x^{\prime}\in\partial S,\ x_{3}\in I\Big{\}},$$

$$\displaystyle\big{(}\mathcal{L}^{h}_{P}(\nabla^{2}_{h}y)\big{)}_{ij}=-{\rm div}_{h}(\partial_{Z}P(\nabla^{2}_{h}y))_{ij*},\ \ \ \ i,j\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\{1,2,3\}$$

$$\displaystyle\begin{cases}-{\rm div}_{h}\Big{(}\partial_{F}W(\nabla_{h}y)+h^{\beta}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{L}^{h}_{P}(\nabla^{2}_{h}y)+\partial_{\dot{F}}R(\nabla_{h}y,\partial_{t}\nabla_{h}y)\Big{)}=h^{3}fe_{3}&\text{in }[0,\infty)\times\Omega\\ y(0,\cdot)=y_{0}&\text{in }\Omega\\ y(t,\cdot)\in\mathfrak{S}_{h}&\text{for }t\in[0,\infty)\end{cases}$$

$$Q_W^2(G)=\underset{a\in {\mathbb{R}}^3}{\min }{Q}_W^3(G^{\ast }+a\otimes e_3+e_3\otimes a)$$

$$Q_D^3(F)=\frac{1}{2}{\partial }_{F_1^2}^2{D}^2(\mathbf{Id},\mathbf{Id})[F,F],Q_D^2(G)=\underset{a\in {\mathbb{R}}^3}{\min }{Q}_D^3(G^{\ast }+a\otimes e_3+e_3\otimes a).$$

$$Q_W^2(G)={\mathbb{C}}_W^2[G,G],Q_D^2(G)={\mathbb{C}}_D^2[G,G]\forall G\in {\mathbb{R}}^{2\times 2}.$$

$$\displaystyle u(x^{\prime}):=\frac{1}{h^{2}}\int_{I}\Big{(}\begin{pmatrix}y_{1}\\ y_{2}\end{pmatrix}(x^{\prime},x_{3})-\begin{pmatrix}x_{1}\\ x_{2}\end{pmatrix}\Big{)}\,dx_{3},\ \ \ \ \ v(x^{\prime}):=\frac{1}{h}\int_{I}y_{3}(x^{\prime},x_{3})\,dx_{3},$$

$$\mathcal{S}=\{(u,v)\in W_0^{1,2}(S;{\mathbb{R}}^2)\times W_0^{2,2}(S)\}.$$

$$\displaystyle\begin{cases}{\rm div}\Big{(}\mathbb{C}^{2}_{W}\big{(}e(u)+\frac{1}{2}\nabla v\otimes\nabla v\big{)}+\mathbb{C}^{2}_{D}\big{(}e(\partial_{t}u)+\nabla\partial_{t}v\odot\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla v\big{)}\Big{)}=0,&\vspace{0.1cm}\textstyle\\ -{\rm div}\Big{(}\Big{(}\mathbb{C}^{2}_{W}\big{(}e(u)+\frac{1}{2}\nabla v\otimes\nabla v\big{)}+\mathbb{C}^{2}_{D}\big{(}e(\partial_{t}u)+\nabla\partial_{t}v\odot\nabla v\big{)}\Big{)}\nabla v\Big{)}&\vspace{0.1cm}\textstyle\\ \quad\quad\quad\quad\quad+\tfrac{1}{12}{\rm div}\,{\rm div}\Big{(}\mathbb{C}^{2}_{W}\nabla^{2}v+\mathbb{C}^{2}_{D}\nabla^{2}\partial_{t}v\Big{)}=f&\text{in }[0,\infty)\times S\\ u(0,\cdot)=u_{0},\ v(0,\cdot)=v_{0}&\text{in }S\\ (u(t,\cdot),v(t,\cdot))\in{\mathscr{S}}_{0}&\text{ for }t\in[0,\infty)\end{cases}$$

$$\displaystyle\int_{S}\Big{(}\mathbb{C}^{2}_{W}\big{(}e(u)+\tfrac{1}{2}\nabla v\otimes\nabla v\big{)}+\mathbb{C}^{2}_{D}\big{(}e(\partial_{t}u)+\nabla\partial_{t}v\odot\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla v\big{)}\Big{)}:\nabla\varphi_{u}=0,$$

$$\displaystyle\int_{S}\Big{(}\mathbb{C}^{2}_{W}\big{(}e(u)+\tfrac{1}{2}\nabla v\otimes\nabla v\big{)}\Big{)}:\big{(}\nabla v\odot\nabla\varphi_{v}\big{)}$$

$$\displaystyle\quad\quad\quad\quad+\int_{S}\Big{(}\mathbb{C}^{2}_{D}\big{(}e(\partial_{t}u)+\nabla\partial_{t}v\odot\nabla v\big{)}\Big{)}:\big{(}\nabla v\odot\nabla\varphi_{v}\big{)}$$

$$\displaystyle\quad\quad\quad\quad+\frac{1}{12}\int_{S}\Big{(}\mathbb{C}^{2}_{W}\nabla^{2}v+\mathbb{C}^{2}_{D}\nabla^{2}\partial_{t}v\Big{)}:\nabla^{2}\varphi_{v}=\int_{S}f\varphi_{v},$$

$$\displaystyle{\phi}(u,v):=\int_{S}\frac{1}{2}Q_{W}^{2}\Big{(}e(u)+\frac{1}{2}\nabla v\otimes\nabla v\Big{)}+\frac{1}{24}Q_{W}^{2}(\nabla^{2}v)-\int_{S}fv$$

$$\mathcal{D}((u_0,v_0),(u_1,v_1))$$

$$\displaystyle\ \ \ \ \ +\frac{1}{12}Q_{D}^{2}\big{(}\nabla^{2}v_{1}-\nabla^{2}v_{0}\big{)}\Big{)}^{1/2}$$

$$\mathcal{D}(u(s),u(t))\le {\int }_s^{t}m(r)dr\text{for all}a\le s\le t\le b.$$

$$|u^{\prime }{|}_{\mathcal{D}}(t):=\underset{s\to t}{\lim }\frac{\mathcal{D}(u(s),u(t))}{|s-t|}.$$

$$|\varphi (v(t))-\varphi (v(s))|\le {\int }_s^{t}g(v(r))|v^{\prime }{|}_{\mathcal{D}}(r)dr\text{for all}a<s\le t<b.$$

$$|\partial \varphi {|}_{\mathcal{D}}(u):=\underset{w\to u}{\mathrm{limsup}}\frac{(\varphi (u)-\varphi (w){)}^{+}}{\mathcal{D}(u,w)}.$$

Full text

∗ Corresponding author: Martin Kružík

Manuel Friedrich

Institute for Computational and Applied Mathematics

University of Münster

Einsteinstr. 62, D-48149 Münster, Germany

Martin Kružík∗

Czech Academy of Sciences, Institute of Information Theory and Automation

Pod vodárenskou věží 4, CZ-182 08 Praha 8, Czechia

Faculty of Civil Engineering, Czech Technical University

Thákurova 7, CZ-166 29 Praha 6, Czechia

Jan Valdman

Czech Academy of Sciences, Institute of Information Theory and Automation

Pod vodárenskou věží 4, CZ-182 08 Praha 8, Czechia

Institute of Mathematics, University of South Bohemia

Branišovská 1760, CZ-370 05 České Budějovice, Czechia

This paper is dedicated to Alexander Mielke in the occasion of his 60th birthday.

1. Introduction

Neglecting inertia, a nonlinear viscoelastic material in Kelvin’s-Voigt’s rheology (i.e., a spring and a dashpot coupled in parallel) obeys the following system of equations

[TABLE]

Here, $[0,T]$ is the process time interval with $T>0$, $\Omega\subset\mathbb{R}^{3}$ is a smooth bounded domain representing the reference configuration, and $y:[0,T]\times\Omega\to\mathbb{R}^{3}$ is the deformation mapping with corresponding deformation gradient $\nabla y$. Further, W:\mathbb{R}^{3\times 3}\to[0,\infty]\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} is a stored energy density, which represents a potential of the first Piola-Kirchhoff stress tensor ${T^{E}}$, i.e., ${T^{E}}:=\partial_{F}W$, and $F\in\mathbb{R}^{3\times 3}$ is the placeholder for $\nabla y$. Moreover, R:\mathbb{R}^{3\times 3}\times\mathbb{R}^{3\times 3}\to[0,\infty)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} denotes a (pseudo)potential of dissipative forces, where $\dot{F}\in\mathbb{R}^{3\times 3}$ is the placeholder of $\partial_{t}\nabla y$. Finally, $f:\Omega\to\mathbb{R}^{3}$ is a volume density of external forces acting on $\Omega$.

A standard assumption for $W$ is frame indifference, i.e., $W(F)=W(QF)$ for every proper rotation $Q\in{\rm SO}(3)$ and every $F\in\mathbb{R}^{3\times 3}$. This implies that $W$ depends on the right Cauchy-Green strain tensor $C:=F^{\top}F$, see e.g. [17]. The second term on the left-hand side of (1) is the stress tensor ${S}(F,\dot{F}):=\partial_{\dot{F}}R(F,\dot{F})$ which has its origin in viscous dissipative mechanisms of the material. Notice that its potential $R$ plays an analogous role as $W$ in the case of purely elastic, i.e., non-dissipative processes. Naturally, we require that $R(F,\dot{F})\geq R(F,0)=0$. The viscous stress tensor must comply with the time-continuous frame-indifference principle, meaning that ${S}(F,\dot{F})=F\tilde{S}(C,\dot{C})$, where $\tilde{S}$ is a symmetric matrix-valued function. This condition constraints $R$ so that [5, 6, 31]

[TABLE]

for some nonnegative function $\tilde{R}$. In other words, $R$ must depend on the right Cauchy-Green strain tensor $C$ and its time derivative $\dot{C}$.

Recently, in [22], the first two authors proved the existence of weak solutions to equations of the form (1) in three-dimensional nonlinear viscoelasticity for nonsimple materials. While the elastic properties of simple elastic materials depend only on the first gradient, the notion of a nonsimple (or second-grade) material refers to the fact that the elastic energy additionally depends on the second gradient of the deformation. This concept, pioneered by Toupin [42, 43], has proved to be useful in modern mathematical elasticity, see e.g. [8, 9, 15, 21, 32, 33, 38]. Adopting this setting currently appears to be inevitable to establish the existence of solutions, see [22], and [31] for a general discussion about the interplay between the elastic energy and viscous dissipation. We emphasize, however, that a main justification of the investigated model is the observation that, in the small strain limit, the problem leads to the standard system of linear viscoelasticity without second gradient.

In the present paper, we are interested in the analysis of lower-dimensional analogs of (1) which are derived by considering (1) for thin viscoelastic plates and by passing to the vanishing-thickness limit. Such studies, often referred to as dimension reduction, play a significant role in nonlinear analysis and numerics since they allow for simpler computational approaches still preserving main features of the full-dimensional system. In particular, it is important that the relationship between the original models and their lower-dimensional counterparts is made rigorous. Usually, the main tools in a variational setting are $\Gamma$-convergence [19] and geometric rigidity estimates [24]. We refer to [29, 30] for a derivation of membrane models from three-dimensional elasticity or to [16, 24, 25, 36] for analogous approaches to plate theory.

In the framework of nonsimple viscoelastic materials, such a scenario was recently studied by the first two authors in [23], where a von Kármán-like viscoelastic plate model has been identified as an effective 2D dimension-reduction limit. For this analysis, besides rigidity estimates and $\Gamma$-convergence, the main tools are gradient flows in metric spaces developed in [3, 34, 40, 41]. Although there are previous works on viscoelastic plates [10, 37], some even including inertial effects [11, 12], their starting point is already a plate model. In contrast, [23] provides a rigorous derivation from a three-dimensional model of viscoelasticity at finite strains by (i) showing the existence of solutions to the effective 2D system, and by (ii) by proving that these solutions are in a certain sense the limits of solutions to the 3D equations for vanishing thickness.

The main aim of this contribution is to carry out a finite-element convergence analysis of a fully discrete viscoelastic plate model and to investigate its behavior by computational experiments. As a byproduct, we also obtain an alternative existence proof for solutions to the effective 2D system. This analysis is based on proving an abstract convergence result of time-space discretizations to metric gradient flows, see Theorem 3.2.

At many spots, our strategy relies on results obtained in [23] and on the theory of gradient flows in metric spaces [3] which provide us with a robust approach to quasistatic evolutionary problems. In particular, Theorem 3.2 exploits a sequence of minimization problems to construct fully discrete approximations (see (19) and (5)) of curves of maximal slope which are then solutions to the viscoelastic plate equations. This makes the proof partially constructive and, at the same time, it suggests a numerical method to be used.

The plan of the paper is as follows. Section 2 reviews equations of nonlinear viscoelasticity in the framework of nonsimple materials and the resulting system for the von Kármán plates. Mathematical tools from the theory of gradient flows in metric spaces [3], such as generalized minimizing movements and curves of maximal slope [2, 20], are introduced in Section 3. Moreover, Section 3 contains our main abstract convergence result for time-space discretizations whose limits are curves of maximal slope, see Theorem 3.2. Section 4 applies the abstract results to the 2D system of viscoelastic von Kármán plate equations, see (15) below: we provide an approximation of the original problem by a finite element method, see Theorem 4.1. As a byproduct, this approximation result yields an alternative proof of the existence of solutions to the viscoelastic plate model originally obtained in [23]. Finally, Section 5 provides computational examples simulating the behavior of the viscoelastic plate exposed to external forces.

We use standard notation for Lebesgue spaces, $L^{p}(\Omega)$, which consist of measurable maps on $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, that are integrable with the $p$-th power (if $1\leq p<+\infty$) or essentially bounded (if $p=+\infty$). With $W^{k,p}(\Omega)$ we denote Sobolev spaces, i.e., linear spaces of maps which, together with their weak derivatives up to the order $k\in\mathbb{N}$, belong to $L^{p}(\Omega)$. Further, $W^{k,p}_{0}(\Omega)$ contains maps from $W^{k,p}(\Omega)$ having zero boundary conditions (in the sense of traces). To emphasize the target space $\mathbb{R}^{k}$, $k=1,2,3$, we write $L^{p}(\Omega;\mathbb{R}^{k})$. If $k=1$, we write $L^{p}(\Omega)$ as usual. We refer to [1] for more details on Sobolev spaces. We also denote the components of vector functions $y$ by $y_{1}$, $y_{2}$, and $y_{3}$, and so on. By $\mathbf{Id}$ we denote the identity matrix in $\mathbb{R}^{3\times 3}$. If $A\in\mathbb{R}^{3\times 3\times 3\times 3}$ and $e\in\mathbb{R}^{3\times 3}$, then $Ae\in\mathbb{R}^{3\times 3}$ is such that for $i,j\in\{1,2,3\}$ we define $(Ae)_{ij}:=A_{ijkl}e_{kl}$ where we use Einstein’s summation convention. An analogous convention is used in similar occasions, in the sequel. Finally, at many spots, we closely follow the notation introduced in [3] to ease readability of our work because the theory developed there is one of the main tools of our analysis.

2. Equations of viscoelasticity in 3D and 2D

We first introduce a 3D setting following the setup in [23, 25]. We consider a right-handed orthonormal system $\{e_{1},e_{2},e_{3}\}$ and $S\subset\mathbb{R}^{2}$ open, bounded with Lipschitz boundary, in the span of $e_{1}$ and $e_{2}$. Let $h>0$ small. We consider deformations $w:S\times(-\frac{h}{2},\frac{h}{2})\to\mathbb{R}^{3}$. It is convenient to work in a fixed domain $\Omega=S\times I$ with $I:=(-\frac{1}{2},\frac{1}{2})$ and to rescale deformations according to $y(x)=w(x^{\prime},hx_{3})$, so that $y:\Omega\to\mathbb{R}^{3}$, where we use the abbreviation $x^{\prime}=(x_{1},x_{2})$. We also introduce the notation $\nabla^{\prime}y=y_{,1}\otimes e_{1}+y_{,2}\otimes e_{2}$ for the in-plane gradient, and the scaled gradient

[TABLE]

Moreover, we define the scaled second gradient by

[TABLE]

where $\delta_{3j},\delta_{3k}$ denotes the Kronecker delta.

Stored elastic energy density and body forces: We assume that $W:\mathbb{R}^{3\times 3}\to[0,\infty]$ is a single-well, frame-indifferent stored energy density with the usual assumptions in nonlinear elasticity. We suppose that there exists $c>0$ such that

[TABLE]

where $SO(3)=\{Q\in\mathbb{R}^{3\times 3}:Q^{\top}Q=\mathbf{Id},\,\det Q=1\}$. Moreover, for $p>3$, let $P:\mathbb{R}^{3\times 3\times 3}\to[0,\infty]$ be a higher order perturbation satisfying

[TABLE]

for $0<c_{1}<c_{2}$. Finally, $f\in L^{\infty}(\Omega)$ denotes a volume normal force, i.e., a force oriented in the $e_{3}$ direction.

Dissipation potential and viscous stress: We now introduce a dissipation potential. We follow here the discussion in [31, Section 2.2] and [23, Section 2]. Consider a time-dependent deformation $y:[0,T]\times\Omega\to\mathbb{R}^{3}$. Viscosity is not only related to the strain rate $\partial_{t}\nabla_{h}y(t,x)$ but also to the strain $\nabla_{h}y(t,x)$. It can be expressed in terms of a dissipation potential $R(\nabla_{h}y,\partial_{t}\nabla_{h}y)$, where $R:\mathbb{R}^{3\times 3}\times\mathbb{R}^{3\times 3}\to[0,\infty)$. An admissible potential has to satisfy frame indifference in the sense (see [5, 31])

[TABLE]

for all $F\in GL_{+}(3)$ and $\dot{F}\in\mathbb{R}^{3\times 3}$, where $GL_{+}(3)=\{F\in\mathbb{R}^{3\times 3}:\det F>0\}$ and $\mathbb{R}^{3\times 3}_{\rm skew}=\{A\in\mathbb{R}^{3\times 3}:A=-A^{\top}\}$.

From the viewpoint of modeling, it is more convenient to postulate the existence of a (smooth) global distance $D:GL_{+}(3)\times GL_{+}(3)\to[0,\infty)$ satisfying $D(F,F)=0$ for all $F\in GL_{+}(3)$. From this, an associated dissipation potential $R$ can be calculated by

[TABLE]

for $F\in GL_{+}(3)$ and $\dot{F}\in\mathbb{R}^{3\times 3}$. Here, $\partial^{2}_{F_{1}^{2}}D^{2}(F_{1},F_{2})$ denotes the Hessian of $D^{2}$ in the direction of $F_{1}$ at $(F_{1},F_{2})$, which is a fourth order tensor. For some $c>0$ we suppose that $D$ satisfies

[TABLE]

Note that conditions (i)-(iii) state that $D$ is a true distance when restricted to symmetric matrices with nonnegative determinants. We cannot expect more due to the separate frame indifference (v). We also point out that (v) implies (6) as shown in [31, Lemma 2.1]. Note that in our model we do not require any conditions of polyconvexity [7] neither for $W$ nor for $D$. One possible example of $D$ satisfying (2) is D(F_{1},F_{2})=|F_{1}^{\top}F_{1}-F_{2}^{\top}F_{2}|\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}. This choice leads to R(F,\dot{F})=|F^{\top}\dot{F}+\dot{F}^{\top}F\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|^{2}/2. For further examples we refer to [31, Section 2.3].

Equations of viscoelasticity in a rescaled domain: Following the study in [23, 28], we introduce the set of admissible configurations by

[TABLE]

where $I=(-\frac{1}{2},\frac{1}{2})$. Note that in [23, 28] more general clamped boundary conditions are considered that are not included here for the sake of simplicity. We formulate the equations of viscoelasticity for a nonsimple material involving the perturbation $P$ (cf. (5)). We introduce a differential operator associated with $P$. To this end, we recall the notation of the scaled gradients in (2)-(3). For i,j\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\{1,2,3\}, we denote by $(\partial_{Z}P(\nabla^{2}_{h}y))_{ij*}$ the vector-valued function $((\partial_{Z}P(\nabla^{2}_{h}y))_{ijk})_{k=1,2,3}$. We also introduce the scaled (distributional) divergence ${\rm div}_{h}g$ for a function $g\in L^{1}(\Omega;\mathbb{R}^{3})$ by ${\rm div}_{h}g=\partial_{1}g_{1}+\partial_{2}g_{2}+\frac{1}{h}\partial_{3}g_{3}$. We define

[TABLE]

for $y\in\mathfrak{S}_{h}$. Let $0<\beta<4$. The equations of nonlinear viscoelasticity are defined by

[TABLE]

for some $y_{0}\in\mathfrak{S}_{h}$, where $\partial_{F}W(\nabla_{h}y)$ $+h^{\beta}\mathcal{L}^{h}_{P}(\nabla^{2}_{h}y)$ denotes the first Piola-Kirchhoff stress tensor and $\partial_{\dot{F}}R(\nabla_{h}y,\partial_{t}\nabla_{h}y)$ the viscous stress with $R$ as introduced in (7).

We remark that the scaling of the forces corresponds to the so-called von Kármán regime. The choice $0<\beta<4$ ensures that the second-gradient term in the energy vanishes in the effective 2D limiting model as $h\to 0$.

Quadratic forms: To formulate the effective 2D problem, we need to consider various quadratic forms. First, we define $Q_{W}^{3}:\mathbb{R}^{3\times 3}\to\mathbb{R}$ by $Q_{W}^{3}(F)=\partial^{2}_{F^{2}}W(\mathbf{Id})[F,F]$. One can show that it depends only on the symmetric part $\frac{1}{2}(F^{\top}+F)$ and that it is positive definite on $\mathbb{R}^{3\times 3}_{\rm sym}=\{A\in\mathbb{R}^{3\times 3}:A=A^{\top}\}$. We also introduce $Q_{W}^{2}:\mathbb{R}^{2\times 2}\to\mathbb{R}$ by

[TABLE]

for $G\in\mathbb{R}^{2\times 2}$, where the entries of $G^{*}\in\mathbb{R}^{3\times 3}$ are given by $G^{*}_{ij}=G_{ij}$ for $i,j\in\{1,2\}$ and zero otherwise. Note that (11) corresponds to a minimization over stretches in the $e_{3}$ direction. In [23] it was assumed that the minimum in (11) is attained for $a=0$. Similarly, we define

[TABLE]

We again assume that the minimum is attained for $a=0$. The assumption that $a=0$ is a minimum in (11)-(12) corresponds to a model with zero Poisson’s ratio in the $e_{3}$ direction. This assumption is not needed in the purely static analysis [25, 28]. However, it is adopted in [23] to simplify the study of the evolutionary problem. We also introduce corresponding symmetric fourth order tensors $\mathbb{C}^{2}_{W}$ and $\mathbb{C}^{2}_{D}$ by

[TABLE]

One can check that $Q^{2}_{W}$ and $Q^{2}_{D}$ are positive semi-definite, and positive definite on $\mathbb{R}_{\rm sym}^{2\times 2}$.

Equations of viscoelasticity in 2D: We now present the effective 2D equations which are formulated in terms of in-plane and out-of-plane displacements fields $u$ and $v$. Following the discussion in [25], these displacement fields can be related to the deformation $y$ in the three-dimensional setting by

[TABLE]

where again $I=(-\frac{1}{2},\frac{1}{2})$. Let us consider the set of admissible displacement fields

[TABLE]

(Compare with (9).) From now on, we are going to work exclusively on the domain $S\subset\mathbb{R}^{2}$ and therefore $\nabla$ will denote the gradient with respect to $x_{1}$ and $x_{2}$, i.e., we will drop the apostrophe from the notation.

Given $(u_{0},v_{0})\in\mathscr{S}$, we consider the equations

[TABLE]

where $\mathbb{C}^{2}_{W}$ and $\mathbb{C}^{2}_{D}$ are defined in (13), and $\odot$ denotes the symmetrized tensor product. Note that the frame indifference of the energy and the dissipation (see (4)(ii) and (2)(v), respectively) imply that the contributions only depend on the symmetric part of the strain $e(u):=\frac{1}{2}(\nabla u+(\nabla u)^{\top})$ and the strain rate $e(\partial_{t}u):=\frac{1}{2}(\partial_{t}\nabla u+\partial_{t}(\nabla u)^{\top})$. Here, ${\rm div}$ denotes the distributional divergence in dimension two.

We also say that $(u,v)\in W^{1,2}([0,\infty);{\mathscr{S}})$ is a weak solution of (15) if $u(0,\cdot)=u_{0}$, $v(0,\cdot)=v_{0}$ and for a.e. $t\geq 0$ we have

[TABLE]

for all $\varphi_{u}\in W^{1,2}_{0}(S;\mathbb{R}^{2})$ and $\varphi_{v}\in W^{2,2}_{0}(S)$. Note that (16a) corresponds to two and (16b) corresponds to one equation, respectively. It is proved in [23, Thm. 2.2 and Thm. 2.3] that solutions to a semidiscretized-in-time system (10) converge to weak solutions (in the sense of (16)) to the initial-boundary value problem (15).

The following von Kármán energy functional $\phi:\mathscr{S}\to\mathbb{R}$ and the global dissipation distance $\mathcal{D}:\mathscr{S}\times\mathscr{S}\to\mathbb{R}$ due to viscosity will play an important role in our analysis: we define

[TABLE]

for $(u,v)\in{\mathscr{S}}$ and

[TABLE]

for $(u_{0},v_{0}),(u_{1},v_{1})\in{\mathscr{S}}$.

In the next sections, we provide mathematical tools which will be used to show that fully discretized solutions to (15) (i.e., discretized in time and space) converge to weak solutions (in the sense of (16)) to (15).

3. An abstract convergence result

In this section we first recall the relevant definitions for metric gradient flows. Then, based on [34], we prove an abstract convergence result of time-space discretizations to curves of maximal slope.

3.1. Definitions: Curves of maximal slope and time-discrete solutions

We consider a complete metric space $(\mathscr{S},\mathcal{D})$. We say a curve $u:(a,b)\to\mathscr{S}$ is absolutely continuous with respect to $\mathcal{D}$ if there exists $m\in L^{1}(a,b)$ such that

[TABLE]

The smallest function $m$ with this property, denoted by $|u^{\prime}|_{\mathcal{D}}$, is called the metric derivative of $u$ and satisfies for a.e. $t\in(a,b)$ (see [3, Theorem 1.1.2] for the existence proof)

[TABLE]

We now define the notion of a curve of maximal slope. We only give the basic definition here and refer to [3, Section 1.2, 1.3] for motivations and more details. By $h^{+}:=\max(h,0)$ we denote the positive part of a function $h$.

Definition 3.1** (Upper gradients, slopes, curves of maximal slope).**

We consider a complete metric space $(\mathscr{S},\mathcal{D})$ with a functional $\phi:\mathscr{S}\to(-\infty,+\infty]$.

(i) A function $g:\mathscr{S}\to[0,\infty]$ is called a strong upper gradient for $\phi$ if for every absolutely continuous curve $v:(a,b)\to\mathscr{S}$ the function $g\circ v$ is Borel and

[TABLE]

(ii) For each $u\in\mathscr{S}$ the local slope of $\phi$ at $u$ is defined by

[TABLE]

(iii) An absolutely continuous curve $u:(a,b)\to\mathscr{S}$ is called a curve of maximal slope for $\phi$ with respect to the strong upper gradient $g$ if for a.e. $t\in(a,b)$

[TABLE]

We introduce time-discrete solutions for a functional $\phi:\mathscr{S}\to(-\infty,+\infty]$ and the metric $\mathcal{D}$ by solving suitable time-incremental minimization problems: consider a fixed time step $\tau>0$ and suppose that an initial datum $Y^{0}_{\tau}$ is given. Whenever $Y_{\tau}^{0},\ldots,Y^{n-1}_{\tau}$ are known, $Y^{n}_{\tau}$ is defined as (if existent)

[TABLE]

We suppose that for a choice of $\tau$ a sequence $(Y_{\tau}^{n})_{n\in\mathbb{N}}$ solving (19) exists. Then we define the piecewise constant interpolation by

[TABLE]

We call $\tilde{Y}_{\tau}$ a time-discrete solution. Note that the existence of such solutions is usually guaranteed by the direct method of the calculus of variations under suitable compactness, coercivity, and lower semicontinuity assumptions.

3.2. Curves of maximal slope as limits of time-space discretizations

In this subsection we formulate a result about the approximation of curves of maximal slope. It is based on a result in [34] recalled in Subsection 3.3 below. We first state our assumptions. We again consider a complete metric space $(\mathscr{S},\mathcal{D})$ and a functional $\phi:\mathscr{S}\to[0,\infty]$. Although $\mathcal{D}$ naturally induces a topology on $\mathscr{S}$, it is often convenient to consider a weaker Hausdorff topology $\sigma$ on $\mathscr{S}$ to have more flexibility in the derivation of compactness properties (see [3, Remark 2.0.5]). We assume that for each $n\in\mathbb{N}$ there exists a $\sigma$-sequentially compact set $K_{N}\subset\mathscr{S}$ such that

[TABLE]

Moreover, we suppose that the topology $\sigma$ satisfies

[TABLE]

We further assume the existence of mutual recovery sequences: for each sequence $z_{k}\stackrel{{\scriptstyle\sigma}}{{\to}}z$ and $w\in\mathscr{S}$ there exists a sequence $(w_{k})_{k}\subset\mathscr{S}$ such that

[TABLE]

This condition is reminiscent of [32, (2.1.37)]. We also point out that this assumption is weaker than the one considered in [34, (2.26)-(2.27)].

We consider a sequence of subspaces $\mathscr{S}_{k}\subset\mathscr{S}$, $k\in\mathbb{N}$, such that each $\mathscr{S}_{k}$ is closed with respect to the topology $\sigma$. By $\rho$ we denote a stronger topology on $\mathscr{S}$ with the property that $\phi$ and $\mathcal{D}$ are continuous with respect to $\rho$. We suppose that $\bigcup_{k}\mathscr{S}_{k}$ is $\rho$-dense in $\mathscr{S}$, i.e., for each $z\in\mathscr{S}$ we find a sequence \color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(z_{k})_{k}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\mathscr{S} such that

[TABLE]

In our applications, $\mathscr{S}_{k}$ will represent finite element subspaces.

Finally, we require a property about geodesical convexity: let $M>0$ and let $\Theta^{1}_{M},\Theta^{2}_{M}:[0,+\infty)\to[0,+\infty)$ be continuous, increasing functions which satisfy $\lim_{t\to 0}\Theta^{1}_{M}(t)/t=1$ and $\lim_{t\to 0}\Theta^{2}_{M}(t)/t=0$. We suppose that for all $z_{0},z_{1}\in\mathscr{S}$ with $\phi(z_{0})\leq M$ there exists a curve $(\gamma_{s})_{s\in[0,1]}\subset\mathscr{S}$ with $\gamma_{0}=z_{0}$ and $\gamma_{1}=z_{1}$ such that

[TABLE]

Moreover, we assume that, if $z_{0},z_{1}$ in $\mathscr{S}_{k}$, then $(\gamma_{s})_{s\in[0,1]}\subset\mathscr{S}_{k}$, as well. In our applications, these curves will simply be convex combinations and, in this context, we will exploit that the finite element spaces $\mathscr{S}_{k}$ are obviously convex sets.

Remark 1** (Convexity assumption on $\Phi$).**

We mention that the condition presented here is tailor-made for our applications to the viscoelastic plate model since in this case we can find curves satisfying (3.2) for specific $\Theta_{M}^{1}$ and $\Theta_{M}^{2}$, see Lemma 4.4 below. We point out that the condition is slightly more general than the one used in [3, Assumption 2.4.5] or [34, Assumption 9]: fix $\lambda\in\mathbb{R}$ and let $\lambda^{-}=-\lambda$ for $\lambda\leq 0$ and $1/\lambda^{-}=+\infty$ else. We suppose that there exists a curve $(\gamma_{s})_{s\in[0,1]}\subset\mathscr{S}$ with $\gamma_{0}=z_{0}$ and $\gamma_{1}=z_{1}$ such that for all $\tau\in(0,1/\lambda^{-})$ and all $s\in[0,1]$ there holds

[TABLE]

Note here that the curve $(\gamma_{s})_{s}$ is chosen independently of $\tau$. A prototypical case is the case of $\lambda$-geodesically convex functionals $\phi$, see [3, Definition 2.4.3]. We briefly check that (26) implies (3.2).

In view of (19), multiplying (26) with $2\tau$ and passing to the limit $\tau\to 0$ we obtain

[TABLE]

i.e., (3.2)(i) holds for $\Theta_{M}^{1}(t)=t$. On the other hand, for $\tau\nearrow 1/\lambda^{-}$, we get

[TABLE]

i.e., (3.2)(ii) holds for $\Theta_{M}^{2}(t)=\frac{\lambda^{-}}{2}t^{2}$. In this sense, (26) can be understood as a special case of (3.2) with $\Theta_{M}^{1}$ being the identity and $\Theta_{M}^{2}$ being quadratic.

We now state or main approximation result. Recall the definition of time-discrete solutions $\tilde{Y}_{\tau}$ in (19)-(20). We say that $\tilde{Y}_{\tau}$ is a time-discrete solution in $\mathscr{S}_{k}$ if $\tilde{Y}_{\tau}(0)\in\mathscr{S}_{k}$ and the minimization problem in (19) is restricted to $\mathscr{S}_{k}$.

Theorem 3.2**.**

Let $(\mathscr{S},\mathcal{D})$ be a complete metric space and let $\phi:\mathscr{S}\to[0,\infty]$. Consider topologies $\tau$ and $\rho$ on $\mathscr{S}$ such that $\mathcal{D}$ and $\phi$ are continuous with respect to $\rho$. Consider $\sigma$-closed subspaces $\mathscr{S}_{k}\subset\mathscr{S}$ and suppose that (21)-(3.2) hold. Consider a null sequence $(\tau_{k})_{k}$. Let $\bar{z}_{0}\in\mathscr{S}$.

Then there exist initial values $(Y^{0}_{k,\tau_{k}})_{k}$ satisfying $Y^{0}_{k,\tau_{k}}\in\mathscr{S}_{k}$ and

[TABLE]

sequences of time-discrete solutions $(\tilde{Y}_{k,\tau_{k}})_{k}$ in $\mathscr{S}_{k}$ starting from $(Y^{0}_{k,\tau_{k}})_{k}$, and a limiting curve $z:[0,+\infty)\to\mathscr{S}$ such that up to a subsequence (not relabeled)

[TABLE]

as $k\to\infty$. The function $z$ is a curve of maximal slope for $\phi$ with respect to $|\partial\phi|_{\mathcal{D}}$.

3.3. Curves of maximal slope as limits of time-discrete solutions

In this subsection we recall a result about the limits of time-discrete solutions obtained by Ortner [34] which is the main ingredient for the proof of Theorem 3.2. We consider a set $\mathscr{S}$ and a sequence of metrics $(\mathcal{D}_{k})_{k}$ on $\mathscr{S}$ as well as a limiting metric $\mathcal{D}$. We again assume that all metric spaces are complete. Moreover, let $(\phi_{k})_{k}$ be a sequence of functionals with $\phi_{k}:\mathscr{S}\to[0,\infty]$.

As before, we consider a Hausdorff topology $\sigma$ on $\mathscr{S}$ which is possibly weaker than the one induced by $\mathcal{D}$. We suppose that the topology $\sigma$ satisfies

[TABLE]

Moreover, assume that for all $n\in\mathbb{N}$ there exists a $\sigma$-sequentially compact set $K_{N}\subset\mathscr{S}$ such that for all $k\in\mathbb{N}$

[TABLE]

Specifically, for a sequence $(z_{k})_{k}$ with $\phi_{k}(z_{k})\leq N$, we find a subsequence (not relabeled) and $z\in\mathscr{S}$ such that $z_{k}\stackrel{{\scriptstyle\sigma}}{{\to}}z$. We suppose lower semicontinuity of the energies and the slopes in the following sense: for all $z\in\mathscr{S}$ and sequences $(z_{k})_{k}$, $z_{k}\in\mathscr{S}_{k}$, we have

[TABLE]

We remark that the condition in [34, (2.10)] is slightly stronger than (30b) since there the condition is required for all sequences and not only on sublevel sets of $\phi_{k}$. The following results remain true under the weaker assumption (30b), cf., e.g., [3, Corollary 2.4.12]. Note that nonnegativity of $\phi_{k}$ and $\phi$ can be generalized to a suitable coerciveness condition, see [3, (2.1.2b)] or [34, (2.5)], which we do not include here for the sake of simplicity. We formulate the main convergence result of time-discrete solutions to curves of maximal slope, proved in [34, Section 2].

Theorem 3.3**.**

Suppose that (28)-(30) hold. Moreover, assume that $|\partial\phi|_{\mathcal{D}}$ is a strong upper gradient for $\phi$. Consider a null sequence $(\tau_{k})_{k}$. Let $(Y^{0}_{k,\tau_{k}})_{k}$ with Y^{0}_{k,\tau_{k}}\in\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathscr{S}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} and $\bar{z}_{0}\in\mathscr{S}$ be initial data satisfying

[TABLE]

Then for each sequence of discrete solutions $(\tilde{Y}_{k,\tau_{k}})_{k}$ for $\phi_{k}$ and $\mathcal{D}_{k}$ starting from $(Y^{0}_{k,\tau_{k}})_{k}$, see (19)-(20), there exists a limiting function $z:[0,+\infty)\to\mathscr{S}$ such that up to a subsequence (not relabeled)

[TABLE]

as $k\to\infty$, and $z$ is a curve of maximal slope for $\phi$ with respect to $|\partial\phi|_{\mathcal{D}}$.

For the proof we refer to [34, Proposition 5, 6]. We comment that this convergence result might seem weak at first glance since in the family of approximations there exists only a subsequence converging to a solution. In practice, however, this often does not cause problems, see [34, Remark 7] for a thorough comment.

3.4. Proof of Theorem 3.2

This subsection is devoted to the proof of Theorem 3.2. Consider the complete metric spaces $(\mathscr{S},\mathcal{D})$ and $(\mathscr{S}_{k},\mathcal{D})$, the functional $\phi:\mathscr{S}\to[0,+\infty]$, and recall assumptions (21)-(3.2). We start with a representation of the local slope defined in Definition 3.1. We also define $\phi_{k}:\mathscr{S}\to[0,+\infty]$ by $\phi_{k}(z)=\phi(z)$ if $z\in\mathscr{S}_{k}$ and $\phi_{k}(z)=+\infty$ else.

Lemma 3.4** (Representation of the local slope).**

Let $M>0$. The local slope for the energy ${\phi}$ in the complete metric space $(\mathscr{S},\mathcal{D})$ admits the representation

[TABLE]

for all $z\in{\mathscr{S}}$ with $\phi(z)\leq M$, where $\Theta^{1}_{M}$ and $\Theta^{2}_{M}$ are the functions from (3.2). The local slope is a strong upper gradient for ${\phi}$. The same representation holds for $\phi_{k}$ in place of $\phi$.

Proof.

We prove the result only for $\phi$. The argument for $\phi_{k}$ is exactly the same which we will explain briefly at the end of the proof. We follow the lines of the proofs of Theorem 2.4.9 and Corollary 2.4.10 in [3], see also [23, Lemma 4.9]. Let $M>0$ and $z\in\mathscr{S}$ with $\phi(z)\leq M$. Recall that $\lim_{t\to 0}\Theta^{1}_{M}(t)/t=1$ and $\lim_{t\to 0}\Theta^{2}_{M}(t)/t=0$. We also recall the definition of the local slope in Definition 3.1 and obtain

[TABLE]

In the second equality we used that $w\to z$ means ${\mathcal{D}}(z,w)\to 0$, and the fact that $\lim_{t\to 0}\Theta_{M}^{1}(t)/t=1$ and $\lim_{t\to 0}\Theta^{2}_{M}(t)/t=0$.

To see the other inequality, we fix $w\neq z$. It is not restrictive to suppose that

[TABLE]

Let $(\gamma_{s})_{s\in[0,1]}$ be the curve given in (3.2) with $\gamma_{0}=z$ and $\gamma_{1}=w$. By (3.2) we obtain

[TABLE]

Since $\mathcal{D}(\gamma_{s},z)\to 0$ as $s\to 0$, see (3.2)(i), we conclude

[TABLE]

The claim now follows by taking the supremum with respect to $w\in\mathscr{S}$.

With this representation of the local slope at hand, one can also show that $|\partial{\phi}|_{{\mathcal{D}}}$ is a strong upper gradient. We refer the reader to [3, Corollary 2.4.10] and [23, Lemma 4.9] for details.

The same argument works for $\phi_{k}$ in place of $\phi$. The only important point to notice is that the curve $(\gamma_{s})_{s\in[0,1]}$ lies in $\mathscr{S}_{k}$ if $z,w\in\mathscr{S}_{k}$, see the line below (3.2), i.e., $\phi_{k}(\gamma_{s})=\phi(\gamma_{s})$ for $s\in[0,1]$. ∎

We are now ready for the proof of Theorem 3.2.

Proof of Theorem 3.2.

Consider a null sequence $(\tau_{k})_{k}$ sequence and $\bar{z}_{0}\in\mathscr{S}$. By (24) and the fact that $\phi$ is continuous with respect to $\rho$ (to recall its definition, see the paragraph preceding (24)) we find a sequence $(Y^{0}_{k,\tau_{k}})_{k}$ satisfying $Y^{0}_{k,\tau_{k}}\in\mathscr{S}_{k}$, Y^{0}_{k,\tau_{k}}\stackrel{{\scriptstyle\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\rho\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}}{{\to}}\bar{z}_{0}, and $\phi(Y^{0}_{k,\tau_{k}})\to\phi(\bar{z}_{0})$. This yields (27) and also (3.3)(ii). Since also $\mathcal{D}$ is continuous with respect to $\rho$, (3.3)(i) holds as well.

Recall the definition $\phi_{k}:\mathscr{S}\to[0,+\infty]$ by $\phi_{k}(z)=\phi(z)$ if $z\in\mathscr{S}_{k}$ and $\phi_{k}(z)=+\infty$ else. We define time-discrete solutions $(\tilde{Y}_{k,\tau_{k}})_{k}$ in the sense of (19)-(20) with respect to $\phi_{k}$ starting from $(Y^{0}_{k,\tau_{k}})_{k}$. Their existence follows from the direct method of the calculus of variations, by using (21), (3.2), and the fact that $\mathscr{S}_{k}$ is closed with respect to $\sigma$. Clearly, these correspond to time-discrete solutions for $\phi$ in $\mathscr{S}_{k}$.

It remains to check that the time-discrete solutions converge to a limiting curve which is a curve of maximal slope for $\phi$ with respect to $|\partial\phi|_{\mathcal{D}}$. Our goal is to apply Theorem 3.3. Since (3.3) has already been verified and $|\partial{\phi}|_{{\mathcal{D}}}$ is a strong upper gradient by Lemma 3.4, it remains to confirm (28)-(30). Set $\mathcal{D}_{k}=\mathcal{D}$ for all $k\in\mathbb{N}$. First, (28) and (29) follow from the fact that $\phi\leq\phi_{k}$, (21), and (3.2)(i). In a similar fashion, (30a) follows from (3.2)(ii). We now show (30b).

Consider a sequence $z_{k}\in\mathscr{S}_{k}$ with $\sup_{k}\phi_{k}(z_{k})\leq M<+\infty$ and $z_{k}\stackrel{{\scriptstyle\sigma}}{{\to}}z$. By (30a) we find also $\phi(z)\leq M$. Let $\varepsilon>0$. By applying Lemma 3.4 we choose $w\in\mathscr{S}$ such that

[TABLE]

Let $(w_{k})_{k}\subset\mathscr{S}$ be a mutual recovery sequence as given by (23). By (24) and the fact that $\mathcal{D}$ and $\phi$ are continuous with respect to the topology $\rho$, we can suppose that $w_{k}\in\mathscr{S}_{k}$ and convergence (23) still holds. By (32), $\phi(z_{k})=\phi_{k}(z_{k})$, $\phi(w_{k})=\phi_{k}(w_{k})$, and the fact that $\Theta^{i}_{M}$ is continuous, increasing for $i=1,2$, we then obtain

[TABLE]

By Lemma 3.4 (for $\phi_{k}$) we then get

[TABLE]

As $\varepsilon$ was arbitrary, we get (30b).

The statement now follows from the abstract convergence result formulated in Theorem 3.3. ∎

4. Finite element approximation of weak solutions to von Kármán viscoelastic plates

In this section we apply Theorem 3.2 to our example of von Kármán viscoelastic plates. Let $\mathscr{S}=W^{1,2}_{0}(S;\mathbb{R}^{2})\times W^{2,2}_{0}(S)$, see (14). We denote the strong convergence in $W^{1,2}(S;\mathbb{R}^{2})\times W^{2,2}(S)$ by $\stackrel{{\scriptstyle\rho}}{{\to}}$. Moreover, we introduce a weak topology $\sigma$ on $\mathscr{S}$: we say that $(u_{k},v_{k})\stackrel{{\scriptstyle\sigma}}{{\to}}(u,v)$ if $u_{k}\rightharpoonup u$ weakly in $W^{1,2}(S;\mathbb{R}^{2})$ and $v_{k}\rightharpoonup v$ weakly in $W^{2,2}(S)$. We let $\mathscr{S}_{k}\subset\mathscr{S}$ be finite dimensional subspaces of finite elements such that $\mathscr{S}_{k}$ is closed with respect to $\sigma$ and $\bigcup_{k}\mathscr{S}_{k}$ is $\rho$-dense in $\mathscr{S}$ in the sense of (24). For an example of such spaces we refer to Section 5 below.

We let $\phi$ and $\mathcal{D}$ as defined in (17) and (2), respectively. For simplicity, we set $f\equiv 0$ since the adaptions for the general case are minor and standard.

We recall that $\tilde{Y}_{k,\tau}$ is called a time-discrete solution in $\mathscr{S}_{k}$ if $\tilde{Y}_{k,\tau}(0)\in\mathscr{S}_{k}$ and the minimization problem in (19) is restricted to $\mathscr{S}_{k}$. Our main result is the following.

Theorem 4.1** (Finite element approximation of weak solutions).**

Consider a null sequence $(\tau_{k})_{k}$ and let $(u_{0},v_{0})\in{\mathscr{S}}$. Then there exist initial values $(U^{0}_{k,\tau_{k}})_{k}$, $(V^{0}_{k,\tau_{k}})_{k}$ satisfying $(U^{0}_{k,\tau_{k}},U^{0}_{k,\tau_{k}})_{k}\in\mathscr{S}_{k}$ and

[TABLE]

sequences of time-discrete solutions $(\tilde{U}_{k,\tau_{k}},\tilde{V}_{k,\tau_{k}})_{k}$ in $\mathscr{S}_{k}$ starting from the initial values $(U^{0}_{k,\tau_{k}},U^{0}_{k,\tau_{k}})_{k}$, and a weak solution $(u,v):[0,\infty)\to{\mathscr{S}}$ to the partial differential equations (15) in the sense of (16) such that up to a subsequence (not relabeled)

[TABLE]

as $k\to\infty$.

Note that this theorem provides us with the strong convergence of time-discrete finite-element approximations to a solution to the original problem.

The result relies on our abstract approximation result stated in Theorem 3.2. In order to apply Theorem 3.2, we need to check the assumptions (21)-(3.2). To this end, we recall some of the results obtained in [23].

Lemma 4.2** (Properties of $({\mathscr{S}},{\mathcal{D}})$ and ${\phi}$).**

We have:

• (i)

$({\mathscr{S}},{\mathcal{D}})$* is a complete metric space.*

• (ii)

Compactness: If $(u_{k},v_{k})_{k}\subset{\mathscr{S}}$ is a sequence with $\sup_{k}{\phi}(u_{k},v_{k})<+\infty$, then $(u_{k},v_{k})_{k}$ is bounded in $W^{1,2}(S;\mathbb{R}^{2})\times W^{2,2}(S)$.

• (iii)

Topologies: The topology induced by ${\mathcal{D}}$ is equivalent to the topology $\rho$.

• (iv)

Continuity: ${\mathcal{D}}((u_{k},v_{k}),(u,v))\to 0$    $\Rightarrow$    $\lim_{k\to\infty}{\phi}(u_{k},v_{k})={\phi}(u,v)$.

Proof.

See [23, Lemma 4.6]. ∎

Theorem 4.3** (Curves of maximal slope and weak solutions).**

For all $(u_{0},v_{0})\in{\mathscr{S}}$, each curve of maximal slope $(u,v):[0,\infty)\to{\mathscr{S}}$ for ${\phi}$ with respect to $|\partial{\phi}|_{{\mathcal{D}}}$ with $(u,v)(0)=(u_{0},v_{0})$ is a weak solution to the partial differential equations (15) in the sense of (16).

Proof.

See [23, Theorem 2.2]. ∎

Lemma 4.4** (Convexity and generalized geodesics).**

Let $M>0$. Then there exist smooth increasing functions $\Theta^{1},\Theta^{2}_{M}:[0,\infty)\to[0,\infty)$ satisfying $\lim_{t\to 0}\Theta^{1}(t)/t=1$ and $\lim_{t\to 0}\Theta^{2}_{M}(t)/t=0$ such that for all $(u_{0},v_{0})\in{\mathscr{S}}$ with $\phi(u_{0},v_{0})\leq M$ and all $(u_{1},v_{1})\in{\mathscr{S}}$ there holds

[TABLE]

where $u_{s}:=(1-s)u_{0}+su_{1}$ and $v_{s}:=(1-s)v_{0}+sv_{1}$, $s\in[0,1]$.

Proof.

See [23, Lemma 4.8]. ∎

Lemma 4.5** (Representation of energy and dissipation).**

Let $\Omega=S\times(-\frac{1}{2},\frac{1}{2})$. For $(u,v)\in\mathscr{S}$ we define for brevity

[TABLE]

Then $\phi$ and $\mathcal{D}$ can be represented as

[TABLE]

Proof.

See [23, Remark 5.4] ∎

We are now in a position to prove Theorem 4.1.

Proof of Theorem 4.1.

First, note that $\phi$ and $\mathcal{D}$ are continuous with respect to $\rho$, see Lemma 4.2(iii),(iv). We now check that (21)-(3.2) hold. First, (21) follows from the choice of $\sigma$, Lemma 4.2(ii), and a compactness argument.

Recall (33). Given a sequence $(u_{k},v_{k})_{k}$ with $(u_{k},v_{k})\stackrel{{\scriptstyle\sigma}}{{\to}}(u,v)$, we observe $G(u_{k},v_{k})\rightharpoonup G(u,v)$ weakly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$ since W^{2,2}\subset\subset\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W^{1,4}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0} in dimension two. Then property (3.2) follows from (4.5) and the fact that $Q^{2}_{W}$ and $Q^{2}_{D}$ are positive semi-definite, see below (13). By the definition of $\mathscr{S}_{k}$ we get (24). To see (3.2), we use Lemma 4.4 and the fact that, if $(u_{0},v_{0}),(u_{1},v_{1})\in\mathscr{S}_{k}$, the convex combinations also lie in $\mathscr{S}_{k}$ due to the convexity of the sets $\mathscr{S}_{k}$.

It remains to prove (23). To this end, consider a sequence $(u_{k},v_{k})$ with $(u_{k},v_{k})\stackrel{{\scriptstyle\sigma}}{{\to}}(u,v)$, and recall that $G(u_{k},v_{k})\rightharpoonup G(u,v)$ weakly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$. Suppose that also $(\bar{u},\bar{v})\in\mathscr{S}$ is given. We define $\bar{u}_{k}=u_{k}+\bar{u}-u$ and $\bar{v}_{k}=v_{k}+\bar{v}-v$. Then by (33) and an elementary expansion we get

[TABLE]

Since $v_{k}\to v$ strongly in $W^{1,4}(S)$, we get that $G(u_{k},v_{k})-G(\bar{u}_{k},\bar{v}_{k})$ converges strongly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$ to

[TABLE]

i.e., $G(u_{k},v_{k})-G(\bar{u}_{k},\bar{v}_{k})$ converges strongly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$ to $G(u,v)-G(\bar{u},\bar{v})$. In view of (4.5)(ii), this implies

[TABLE]

Moreover, by an elementary expansion and (4.5)(i) we get

[TABLE]

Since $G(u_{k},v_{k})-G(\bar{u}_{k},\bar{v}_{k})$ converges strongly to $G(u,v)-G(\bar{u},\bar{v})$ and $G(u_{k},v_{k})\rightharpoonup G(u,v)$ weakly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$, we get ${\phi}(u_{k},v_{k})-{\phi}(\bar{u}_{k},\bar{v}_{k})\to{\phi}(u,v)-{\phi}(\bar{u},\bar{v})$. This along with (35) shows that (3.2) holds.

Having checked (21)-(3.2), we can now apply Theorem 3.2. This yields the existence of time-discrete solutions and of a curve of maximal slope such that convergence of time-discrete solutions holds with respect to $\sigma$. The fact that the curve is a weak solution to (15) follows from Theorem 4.3. It remains to prove that the convergence of time-discrete solutions holds with respect to the strong topology $\rho$.

To confirm the latter property, we use the principle that weak convergence together with energy convergence induces strong convergence. More specifically, given $(u_{k},v_{k})\stackrel{{\scriptstyle\sigma}}{{\to}}(u,v)$ and $\phi(u_{k},v_{k})\to\phi(u,v)$, we argue as follows. Since $G(u_{k},v_{k})\rightharpoonup G(u,v)$ weakly in $L^{2}(\Omega;\mathbb{R}^{2\times 2})$, we get by (4.5)(i) that

[TABLE]

Since $Q^{2}_{W}$ is positive definite on $\mathbb{R}_{\rm sym}^{2\times 2}$, see below (13), we get ${\rm sym}(G(u_{k},v_{k}))\to{\rm sym}(G(u,v))$ strongly in $L^{2}(\Omega;\mathbb{R}^{2\times 2}_{\rm sym})$. Then, in view of (33), using Poincaré’s and Korn’s inequality, together with zero boundary conditions, it is elementary to check that $u_{k}\to u$ in $W^{1,2}(S;\mathbb{R}^{2})$ and $v_{k}\to v$ in $W^{2,2}(S)$, i.e., $(u_{k},v_{k})\stackrel{{\scriptstyle\rho}}{{\to}}(u,v)$. ∎

5. Numerical experiments

In this section we describe two numerical experiments on a homogeneous and isotropic viscoelastic plate. Our computational strategy relies on a sequence of minimization problems based on (17) and (2). Take a time horizon $T>0$ and a time step $\tau>0$ such that $n_{max}:=T/\tau\in\mathbb{N}$. Having an initial condition $(u_{0},v_{0})\in\mathscr{S}_{k}$ in a finite element space $\mathscr{S}_{k}$ detailed below, we find, for $1\leq n\leq n_{max}$, a solution $(u_{n},v_{n})\in\mathscr{S}_{k}$ of the following problem

[TABLE]

where $\phi$ and $\mathcal{D}$ are defined in (17) and (2), respectively. As $Q^{2}_{W}$ we take an isotropic material, i.e.,

[TABLE]

and

[TABLE]

for every symmetric $G\in\mathbb{R}^{2\times 2}_{\rm sym}$. The constants $\lambda,\mu$ are Lamé constants and $c>0$ represents a viscosity parameter.

A finite element method (FEM) is applied for the numerical approximation of the in-plane and out-of-plane displacement fields $u$ and $v$. This space corresponds to $\mathscr{S}_{k}$, for some $k\in\mathbb{N}$ large enough, as considered in Section 4. We assume a uniform rectangular mesh in 2D discretizing a square domain

[TABLE]

into square elements with the edge of the length $h=1/k$ and further approximate:

1. (1)

a vector function $u=(u_{1},u_{2})$ by $\text{Q1}_{k}$ elements (elementwise bilinear and globally continuous) in each component $u_{1}$ and $u_{2}$, 2. (2)

a scalar function $v$ by the Bogner-Fox-Schmit ($\text{BFS}_{k}$) rectangular elements [13], i.e., a bi-cubic Hermite elements, that provide globally $C^{1}$ approximations.

We define $\mathscr{S}_{k}=\{(u,v)\in\text{Q1}_{k}\times\text{BFS}_{k}\}$. The $\rho$-denseness of $\bigcup_{k}\mathscr{S}_{k}$ in $\mathscr{S}$ follows from the properties of the finite-element interpolants, see [18, Thms. 3.2.3 and 6.1.7] or [14]. A rectangular mesh with 81 nodes and 64 rectangles is given in Figure 1 together with Gauss integration points used in quadrature formulas.

5.1. Benchmark I

We consider a time sequence of minimization problems with the time step $\tau=1$, Lamé parameters $\lambda=\mu=1e3$, the viscosity parameter $c=3e3$, and the constant volume force $f=-1e3$. The initial condition is given by

[TABLE]

and the boundary condition for $t\geq 0$ by

[TABLE]

Although the choice of parameters above is not physically relevant, the meaning of this benchmark is clear: The initial condition $(u_{0},v_{0})$ is not equilibrated and therefore after some time (the speed of this transition is driven by the viscosity constant $c$) the energy stabilizes at its equilibrium given by the volume force oriented in the gravity direction. The first 8 minimizers are displayed in Figure 2. The scalar field $v$ is displayed as a vertical plate deformation, whereas the vector displacement field $u$ is displayed as a deformed mesh.

5.2. Benchmark II

The theoretical part of this paper covers boundary conditions defined on the full boundary of $S$ only. The computer simulations, however, are possible also for boundary conditions given on a part of the boundary $\partial S$. We consider a time sequence of minimization problems with the time step $\tau=1$, Lamé parameters $\lambda=\mu=1e3$, the viscosity parameter $c=3e3$, and the constant volume force $f=1e2$. The initial condition is given by

[TABLE]

and the boundary condition for $t\geq 0$ by

[TABLE]

so boundary conditions are given on the lower and upper parts of the domain boundary only. The first 8 minimizers are displayed in Figure 3.

5.3. Implementation details

Our Matlab implementation is based on former vectorized codes of [4, 26, 39] that allow for a fast assembly of various finite element matrices. The code is available at

https://www.mathworks.com/matlabcentral/fileexchange/72991

for

download. It includes an own implementation of the Bogner-Fox-Schmit (BFS) rectangular elements for a uniformly refined rectangular mesh, where all rectangular elements are for simplicity of the same size $hx_{1}\times hx_{2}$ (in our computations $hx_{1}=hx_{2}=h$). The basis functions on each rectangle are based on bicubic polynomials, i.e., tensor products of 4 cubic (Hermite) polynomials. They have 16 degrees of freedom with 4 degrees in each of its 4 corner nodes approximating: a function value, its gradient (two components), and the second mixed derivative. Therefore, the initial function $v_{0}$ must have all these fields available in our simulations. Figure 4 depicts $v_{0}$ from (37) represented in terms of BFS elements. We recall that BFS elements were also successfully tested in [27] and their implementation was explained in detail in [35].

Acknowledgments

The research of MF was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044 -390685587, Mathematics Münster: “Dynamics–Geometry–Structure”. MK and JV acknowledge the support by GAČR project 17-04301S.