Continuous/Discontinuous Finite Element Modelling of Kirchhoff Plate Structures in $\mathbb{R}^3$ Using Tangential Differential Calculus

TL;DR

This paper develops a finite element modeling approach for Kirchhoff plates in three-dimensional space using surface differential calculus, incorporating both membrane and bending deformations with a mixed continuity scheme.

Contribution

It introduces a novel finite element formulation employing tangential differential calculus for Kirchhoff plates, including in-plane membrane effects and a mixed continuous/discontinuous element approach.

Findings

Effective modeling of Kirchhoff plates with membrane effects

Finite element method with continuous displacements and discontinuous rotations

Use of $C^0$-elements for plate and membrane discretization

Abstract

We employ surface differential calculus to derive models for Kirchhoff plates including in-plane membrane deformations. We also extend our formulation to structures of plates. For solving the resulting set of partial differential equations, we employ a finite element method based on elements that are continuous for the displacements and discontinuous for the rotations, using $C^0$-elements for the discretisation of the plate as well as for the membrane deformations. Key to the formulation of the method is a convenient definition of jumps and averages of forms that are $d$-linear in terms of the element edge normals.

Figures

Equations

$${\Omega }_t=\{\mathbf{x}\in {\mathbb{R}}^3:|\zeta (\mathbf{x})|<t/2\}$$

$$\mathbf{p}(\mathbf{x})=\mathbf{x}-\zeta (\mathbf{x})\mathbf{n}(\mathbf{x})$$

$$\nabla \mathbf{p}=\mathbf{I}-\zeta \nabla \otimes \mathbf{n}-\mathbf{n}\otimes \mathbf{n}$$

$${\mathbf{P}}_{\Gamma }:=\mathbf{I}-\mathbf{n}\otimes \mathbf{n}$$

$${\nabla }_{\Gamma }:={\mathbf{P}}_{\Gamma }\nabla$$

$$\nabla_{\Gamma}=:\left[\begin{array}[]{>{\displaystyle}c}\frac{\partial}{\partial x_{\Gamma}}\\[8.53581pt] \frac{\partial}{\partial y_{\Gamma}}\\[8.53581pt] \frac{\partial}{\partial z_{\Gamma}}\end{array}\right]$$

$$\bm{v}\otimes\nabla_{\Gamma}:=\left[\begin{array}[]{>{\displaystyle}c>{\displaystyle}c>{\displaystyle}c}\frac{\partial v_{1}}{\partial x_{\Gamma}}&\frac{\partial v_{1}}{\partial y_{\Gamma}}&\frac{\partial v_{1}}{\partial z_{\Gamma}}\\[8.53581pt] \frac{\partial v_{2}}{\partial x_{\Gamma}}&\frac{\partial v_{2}}{\partial y_{\Gamma}}&\frac{\partial v_{2}}{\partial z_{\Gamma}}\\[8.53581pt] \frac{\partial v_{3}}{\partial x_{\Gamma}}&\frac{\partial v_{3}}{\partial y_{\Gamma}}&\frac{\partial v_{3}}{\partial z_{\Gamma}}\end{array}\right]$$

$$\mathbf{u}(\mathbf{x})={\mathbf{u}}_0(\mathbf{p}(\mathbf{x}))-\zeta (\mathbf{x})\mathbf{w}(\mathbf{p}(\mathbf{x}))$$

$$\mathbf{u}=u_n\mathbf{n}-\zeta {\mathbf{P}}_{\Gamma }\mathbf{\theta }$$

$$\mathbf{\epsilon }(\mathbf{\theta }):=\frac{1}{2}\left(\mathbf{\theta }\otimes \nabla +(\mathbf{\theta }\otimes \nabla {)}^{\mathrm{T}}\right)$$

$${\mathbf{e}}_{\Gamma }(\mathbf{\theta }):=\frac{1}{2}\left(\mathbf{\theta }\otimes {\nabla }_{\Gamma }+(\mathbf{\theta }\otimes {\nabla }_{\Gamma }{)}^{\mathrm{T}}\right)$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{\theta })$$

$$={\mathbf{e}}_{\Gamma }(\mathbf{\theta })-({\mathbf{e}}_{\Gamma }(\mathbf{\theta })\cdot \mathbf{n})\otimes \mathbf{n}-\mathbf{n}\otimes ({\mathbf{e}}_{\Gamma }(\mathbf{\theta })\cdot \mathbf{n})$$

$$\mathbf{\theta }={\mathbf{P}}_{\Gamma }\mathbf{\theta }+(\mathbf{\theta }\cdot \mathbf{n})\mathbf{n}$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{\theta })={\mathbf{\epsilon }}_{\Gamma }({\mathbf{P}}_{\Gamma }\mathbf{\theta })+(\mathbf{\theta }\cdot \mathbf{n})\mathbf{\kappa }$$

$$\mathbf{\kappa }:=\nabla \otimes \mathbf{n}$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{\theta })={\mathbf{\epsilon }}_{\Gamma }({\mathbf{P}}_{\Gamma }\mathbf{\theta })$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{u})={\mathbf{\epsilon }}_{\Gamma }(u_n\mathbf{n})-\zeta {\mathbf{\epsilon }}_{\Gamma }({\mathbf{P}}_{\Gamma }\mathbf{\theta })$$

$$\mathbf{\gamma }=\frac{1}{2}\left(2{\mathbf{e}}_{\Gamma }(u_n\mathbf{n})\cdot \mathbf{n}-{\mathbf{P}}_{\Gamma }\mathbf{\theta }\right)$$

$${\mathbf{e}}_{\Gamma }(u_n\mathbf{n})=u_n{\mathbf{e}}_{\Gamma }(\mathbf{n})+\frac{1}{2}\left(\mathbf{n}\otimes {\nabla }_{\Gamma }{u}_n+(\mathbf{n}\otimes {\nabla }_{\Gamma }{u}_n{)}^{\mathrm{T}}\right)$$

$$2{\mathbf{e}}_{\Gamma }(u_n\mathbf{n})\cdot \mathbf{n}={\nabla }_{\Gamma }{u}_n+2{\mathbf{e}}_{\Gamma }(\mathbf{n})\cdot \mathbf{n}{u}_n={\nabla }_{\Gamma }{u}_n$$

$$\mathbf{\gamma }=\frac{1}{2}\left({\nabla }_{\Gamma }{u}_n-{\mathbf{P}}_{\Gamma }\mathbf{\theta }\right)$$

$$\mathbf{u}:=u_n\mathbf{n}-\zeta {\nabla }_{\Gamma }{u}_n$$

$${\mathbf{\epsilon }}_{\Gamma }(u_n\mathbf{n})$$

$$\displaystyle=\frac{1}{2}\bm{P}_{\Gamma}\Big{(}(\nabla_{\Gamma}u_{n})\otimes\bm{n}+\bm{n}\otimes(\nabla_{\Gamma}u_{n})\Big{)}\bm{P}_{\Gamma}$$

$$=0$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{u})=-\zeta {\mathbf{\epsilon }}_{\Gamma }({\nabla }_{\Gamma }{u}_n)$$

$$\mathbf{u}=-\zeta {\nabla }_{\Gamma }{u}_n$$

$$\mathbf{\sigma }=2\mu \mathbf{\epsilon }+\lambda \text{tr}\mathbf{\epsilon }\mathbf{I}$$

$$\mu =\frac{E}{2(1+\nu )},\lambda =\frac{E\nu }{1-{\nu }^2}$$

$${\mathbf{\sigma }}_{\Gamma }:=2\mu {\mathbf{\epsilon }}_{\Gamma }+\lambda \text{tr}{\mathbf{\epsilon }}_{\Gamma }{\mathbf{P}}_{\Gamma }=2\mu {\mathbf{\epsilon }}_{\Gamma }+\lambda {\nabla }_{\Gamma }\cdot \mathbf{u}{\mathbf{P}}_{\Gamma }$$

$${\mathcal{E}}_{\text{P}}$$

Full text

1 Introduction

The Kirchhoff plate model is a fourth order partial differential equation which requires $C^{1}$–continuous elements for constructing conforming finite element methods. To avoid this requirement, nonconforming finite elements can be used; one classical example being the Morley triangle Mo68 which has displacement degrees of freedom in the corner nodes and rotation degrees of freedom at the midpoint of the edges. If we want to solve also for the membrane displacements, it is more straightforward to be able to use only displacement degrees of freedom for both the normal (plate) and tangential (membrane) displacements. To reach this goal, one can instead use the discontinuous Galerkin (dG) method HaLa02 , more efficiently implemented as a $C^{0}$–continuous Galerkin method allowing for discontinuous approximation of derivatives, referred to as the continuous/discontinuous Galerkin, or c/dG, method, first suggested by Engel et al. engel02 , and further developed for plate models by Hansbo et al. HaHeLa10 ; HaHeLa11 ; HaLa03 ; HaLa11 and by Wells and Dung WeDu07 . See also Larsson and Larson LaLa17 for error estimates in the case of the biharmonic problem on a surface. To obtain a continuous model, we combine the plate equation for the normal displacements with the tangential differential equation for the membrane from Hansbo and Larson HaLa14a to obtain a structure with both bending resistance and membrane action. This model is then discretised using continuous finite elements for the membrane and c/dG for the plate, using the same order polynomial in both cases.

The standard engineering approach to constructing plate elements arbitrarily oriented in ${\mathbb{R}}^{3}$ is to use rotation matrices to transform the displacements from a planar element to the actual, common, coordinates, thus transforming the stiffness matrices. In this paper we instead extend the c/dG method to the case of arbitrarily oriented plates, allowing for membrane deformations, directly using Cartesian coordinates in ${\mathbb{R}}^{3}$. We argue that this makes it simpler to implement discrete schemes in general, and in particular the discontinuous Galerkin terms on the element borders. It also gives an analytical model directly expressed in equilibrium equations in physical coordinates.

A particular feature of our method is the handling of the trace terms in the c/dG method. In the recent paper on dG for elliptic problems on smooth surfaces by Dedner, Madhavan, and Stinner DeMaSt13 the definition of the normal to the element faces (tangential to the surface), the conormal, was discussed and different variants tested numerically. In our case, where the surface is piecewise smooth (planar), the definition of the conormal at plate junctures is crucial to the equilibrium. It turns out the proper way to define the jumps and averages of trace quantities that are $d$-linear in the conormal is to compute the trace on the left and right side with the respective unit conormals and adjust the sign on one of the sides with $(-1)^{d}$. This leads to a generalization of the standard jump and averages in the flat case where a fixed conormal is used for both the left and right side in the definition of the jump. Furthermore, the standard formula, where the jump in a product of two functions is represented as the sum of the two products of the averages and jumps of the two functions, also generalizes to this situation. With these tools at hand we may directly use standard discontinuous Galerkin techniques to derive a finite element method for a plate structure. The resulting method takes the same form as a standard c/dG method for a plate. The only difference is the proper definition of jumps and averages. See also JoLaLa17 , where a similar approach was used for the Laplace-Beltrami operator on a surface with sharp edges.

The outline of the paper is as follows: In Section 2 we derive a variational formulation for a plate with arbitrary orientation in $\mathbb{R}^{3}$, in Section 3 we define the relevant traces, including forces and moments, define the averages and jumps of $d$-linear forms, and formulate the interface conditions for a plate structure, in Section 4 we formulate the finite element method, in Section 5 we present numerical examples, and finally we conclude with some remarks in Section 6.

2 Single Plate

2.1 Tangential Differential Calculus

Let $\Gamma$ be a piecewise planar two-dimensional surface imbedded in ${\mathbb{R}}^{3}$, with piecewise constant unit normal $\bm{n}$ and boundary $\partial\Gamma$, split into a Neumann part $\partial\Gamma_{\text{N}}$ where forces and moments are known, and a Dirichlet part $\partial\Gamma_{\text{D}}$ where rotations and displacements are known. For ease of presentation we shall assume that $\partial\Gamma_{\text{N}}=\emptyset$ and that we have zero displacements and rotations on the boundary. The case of $\partial\Gamma_{\text{N}}\neq\emptyset$ is straightforward to implement and will be used in the numerical examples. Mixed boundary conditions are handled equally straightforward.

If we denote the (piecewise) signed distance function relative to $\Gamma$ by $\zeta(\bm{x})$, for $\bm{x}\in\mathbb{R}^{3}$, fulfilling $\nabla\zeta=\bm{n}$, we can define the domain occupied by the shell by

[TABLE]

where $t$ is the thickness of the shell, which for simplicity will be assumed constant. The closest point projection $\bm{p}:\Omega_{t}\rightarrow\Gamma$ is given by

[TABLE]

the Jacobian matrix of which is

[TABLE]

where $\bm{I}$ is the identity and $\otimes$ denotes exterior product. The corresponding linear projector $\bm{P}_{\Gamma}=\bm{P}_{\Gamma}(\bm{x})$, onto the tangent plane of $\Gamma$ at $\bm{x}\in\Gamma$, is given by

[TABLE]

and we can then define the surface gradient $\nabla_{\Gamma}$ as

[TABLE]

The surface gradient thus has three components, which we shall denote by

[TABLE]

For a vector valued function $\bm{v}(\bm{x})$, we define the tangential Jacobian matrix as

[TABLE]

and the surface divergence $\nabla_{\Gamma}\cdot\bm{v}:=\text{tr}\,\bm{v}\otimes\nabla_{\Gamma}$.

2.2 Displacement and Strain

Upon loading, each point $\bm{x}\in\Omega_{t}$, in the plate undergoes a displacement

[TABLE]

where $\bm{u}_{0}$ and $\bm{w}$ are vector fields defined on $\Gamma$, $\bm{u}_{0}$ arbitrary and $\bm{w}$ a tangential vector, $\bm{w}\cdot\bm{n}=0$ on $\Gamma$, or $\bm{w}=\bm{P}_{\Gamma}\bm{\theta}$ with $\bm{\theta}$ arbitrary. Thus, neglecting in-plane extensions for the moment, we can write

[TABLE]

in $\Omega_{t}$. Here $u_{n}=\bm{u}\cdot\bm{n}$.

We introduce the strain tensor $\bm{\varepsilon}$ as

[TABLE]

and define the symmetric part of the tangential Jacobian as

[TABLE]

The in-plane strain tensor $\bm{\varepsilon}_{\Gamma}$ is implemented using the following identity

[TABLE]

If we write

[TABLE]

then

[TABLE]

where

[TABLE]

is the curvature tensor, cf. DeZo95 ; DeZo96 . For planar $\Gamma$, $\bm{n}$ is constant, and this simplifies to

[TABLE]

The total in-plane strain tensor is thus given by

[TABLE]

In DeZo95 ; DeZo96 it is also shown that the mid-plane rotation in the absence of shear deformation is given by $2\bm{e}_{\Gamma}(u_{n}\bm{n})\cdot\bm{n}$, and for shear deformable inextensible shells we thus have the shear deformation vector

[TABLE]

It is is also easy to verify that

[TABLE]

so that, since $\bm{n}\cdot\nabla_{\Gamma}u_{n}=0$,

[TABLE]

since $\bm{n}$ is constant; thus

[TABLE]

In the tangential setting, the Kirchhoff assumption of zero shear deformations can therefore be written

[TABLE]

Furthermore, we find that

[TABLE]

and for inextensible plates we get

[TABLE]

and in this case we thus only obtain contributions to the strain energy from the displacement field

[TABLE]

2.3 Variational Formulations

We shall assume isotropic stress–strain relations,

[TABLE]

where $\bm{\sigma}$ is the stress tensor, and plane stress conditions, for which the Lamé parameters $\lambda$ and $\mu$ are related to Young’s modulus $E$ and Poisson’s ratio $\nu$ via

[TABLE]

For the in-plane stress tensor we find, by projecting (2.29) from left and right,

[TABLE]

The potential energy of the plate is postulated as

[TABLE]

where $\bm{\sigma}:\bm{\varepsilon}=\sum_{ij}\sigma_{ij}\varepsilon_{ij}$ for second order Cartesian tensors $\bm{\sigma}$ and $\bm{\varepsilon}$. Integrating in $\zeta$, we obtain

[TABLE]

Under the assumption of clamped boundary conditions, the corresponding variational problem is to find u_{n}\in H_{0}^{2}(\Gamma)=\{\text{v\in H^{2}(\Gamma)$:$v=\bm{\nu}\cdot\nabla_{\Gamma}v=0$on$\partial\Gamma}\} such that

[TABLE]

for all $v\in H^{2}_{0}(\Gamma)$.

Introducing also membrane deformations, the total potential energy ${\cal E}_{\text{tot}}$ of the plate must take into account both the bending energy ${\cal E}_{\text{P}}$ and the membrane energy ${\cal E}_{\text{M}}$, so that ${\cal E}_{\text{tot}}={\cal E}_{\text{P}}+{\cal E}_{\text{M}}$, where

[TABLE]

Since we wish to use a 3D Cartesian vector field we redefine $\bm{u}:=\bm{u}_{0}$ and $u_{n}:=\bm{n}\cdot\bm{u}$, make use of (2.17), and introduce the function space

[TABLE]

We are then led to the variational problem of finding $\bm{u}\in V$ such that

[TABLE]

for all $\bm{v}\in V$. Introducing the notation

[TABLE]

we may write (2.37) in the more compact form

[TABLE]

For implementation purposes we note that for $\bm{n}$ constant

[TABLE]

and

[TABLE]

2.4 Strong Form

The corresponding strong form of the problem is to find $\bm{u}=u_{n}\bm{n}+\bm{P}_{\Gamma}\bm{u}$ such that

[TABLE]

and

[TABLE]

3 Plate Structures

3.1 Forces and Moments

Consider first a subdomain polygonal subdomain $\omega\subset\Gamma$ of the plate $\Gamma$ with boundary $\partial\omega$ consisting of line segments $\gamma_{i}$. Using Greens formula on $\omega$ we obtain

[TABLE]

where we used the identity $v_{n}=\bm{v}\cdot\bm{n}$ and moved the normal to the first slot in the bilinear form. Letting $\bm{\tau}$ be a unit tangent vector to $\partial\omega$, we may split the last term on the right hand side of (3.2) in normal and tangent contributions as follows

[TABLE]

where the first term is the bending moment. For the second term on the right hand side (3.3), integrating by parts along one of the line segments $\gamma_{i}$, with unit tangent and normal $\bm{\tau}_{i}=\bm{\tau}|_{\gamma_{i}}$ and $\bm{n}_{i}=\bm{n}|_{\gamma_{i}}$, we obtain

[TABLE]

where $\partial\gamma_{i}$ consists of the two end points of the line segment $\gamma_{i}$. We introduce the following notation

[TABLE]

for the normal and tangent components of the force and the moment at each of the line segments $\gamma$ on $\partial\omega$. Furthermore, we introduce the corner, or Kirchhoff, forces

[TABLE]

at a corner $\bm{x}$ associated with a line segment $\gamma_{i}$, which has $\bm{x}$ as one of its endpoints and $\bm{\tau}_{i}$ is the unit tangent vector to $\gamma_{i}$ directed into $\bm{x}$. We then have the identity

[TABLE]

where $\mathcal{X}(\partial\omega)$ is the set of corners on the polygonal boundary $\partial\omega$ and $\mathcal{I}(\bm{x})$ is an enumeration of the two linesegments that has $\bm{x}$ as one of its endpoints.

3.2 Jumps and Averages

Consider a line segment $\gamma$ shared by two plates $\Gamma^{+}$ and $\Gamma^{-}$. We note that the force $\bm{F}^{\pm}$ is an $\mathbb{R}^{3}$ valued 1-form in $\bm{\nu}^{\pm}$ and the moment $M^{\pm}$ is an $\mathbb{R}$ valued 2-form in $\bm{\nu}^{\pm}$. More generally let $w^{\pm}=w^{\pm}(\bm{\nu}^{\pm},\dots\bm{\nu}^{\pm})$ be an $\mathbb{R}^{n}$ valued $d$-linear form in $\bm{\nu}^{\pm}$. Then we define the jump and average at $\gamma$ by

[TABLE]

Note that when both plates $\Gamma^{+}$ and $\Gamma^{-}$ reside in the same plane $\bm{\nu}^{-}=-\bm{\nu}^{+}$ and we recover, using linearity and the simplified notation $w^{\pm}(\bm{\nu}^{\pm},\dots,\bm{\nu}^{\pm})=w(\bm{\nu}^{\pm})$, the standard jump

[TABLE]

and similarly for the average. Finally, let $w^{\pm}_{i}$ be an $\mathbb{R}^{n}$ valued $d_{i}$-linear form in $\bm{\nu}^{\pm}$, then we note that $(w_{1}\cdot w_{2})^{\pm}=w_{1}^{\pm}\cdot w_{2}^{\pm}$ is an $\mathbb{R}$ valued $(d_{1}+d_{2})$-linear form in $\bm{\nu}^{\pm}$ and we have the identity

[TABLE]

where for $n=1$ the scalar product is just usual multiplication of scalars. We may verify (3.16) by

[TABLE]

3.3 Interface Conditions

Consider now a plate structure consisting of a finite number of plates such that at most two plates intersect in a common line segment. For simplicity we consider clamped boundary conditions on the boundary of the structure and focus our attention on the interface conditions at the intersections between the plates. For each line segment $\gamma$ where two plates $\Gamma^{+}$ and $\Gamma^{-}$ intersect we have the interface conditions

[TABLE]

corresponding to continuity of displacements, continuity of the rotation angle, equilibrium of forces, and equilibrium of moments.

Furthermore, at each corner $\bm{x}$, not residing on the boundary of the structure, we require equilibrium of the Kirchhoff forces

[TABLE]

where $\mathcal{I}(\bm{x})$ is an enumeration of the line segments that meet in the corner $\bm{x}$ and $\bm{F}^{\pm}_{\bm{x},i}$ is the Kirchhoff force emanating from plate $\Gamma_{i}^{\pm}$, the two plates that meet in line segment $i$. In other words, there are two contributions associated with each line segment, one for each of the two plates that share the line segment.

4 Finite Element Formulation

4.1 The Mesh and Finite Element Space

Let $\widehat{K}\subset\mathbb{R}^{2}$ be a reference triangle and let $P_{2}(\widehat{K})$ be the space of polynomials of order less or equal to $2$ defined on $\widehat{K}$. Let $\Gamma$ be triangulated with quasi uniform triangulation $\mathcal{K}_{h}$ and mesh parameter $h\in(0,h_{0}]$ such that each triangle $K=F_{K}(\widehat{K})$ is planar (a subparametric formulation). We let $\mathcal{E}_{h}$ denote the set of edges in the triangulation.

We here extend the discontinuous Galerkin method of Dedner et al. DeMaSt13 for the Laplace–Beltrami operator to the case of the plate. We recall that $\Gamma$ is piecewise planar and thus $\bm{n}$ is a piecewise constant exterior unit normal to $\Gamma$.

For the parametrization of $\Gamma$ we wish to define a map from a reference triangle $\widehat{K}$ defined in a local coordinate system $(\xi,\eta)$ to any given triangle $K$ on $\Gamma$. Thus the coordinates of the discrete surface are functions of the reference coordinates inside each element, $\bm{x}_{\Gamma}=\bm{x}_{\Gamma}(\xi,\eta)$. For any given parametrization, we can extend it to $\Omega_{t}$ by defining

[TABLE]

where $-t/2\leq\zeta\leq t/2$ and $\bm{n}$ is the normal to $\Gamma$.

We consider in particular a finite element parametrization of $\Gamma$ as

[TABLE]

where $\bm{x}_{i}$ are the physical location of the (geometry representing) nodes on the initial midsurface and $\psi_{i}(\xi,\eta)$ are affine finite element shape functions on the reference element. (This parametrization is of course exact in the case of a piecewise planar $\Gamma$.)

For the approximation of the displacement, we use a constant extension,

[TABLE]

where $\bm{u}_{i}$ are the nodal displacements, and $\varphi_{i}$ are piecewise quadratic shape functions. We employ the usual finite element approximation of the physical derivatives of the chosen basis $\{\varphi_{i}\}$ on the surface, at $(\xi,\eta)$, in matrix representation, as

[TABLE]

where ${\bm{J}}(\xi,\eta,\zeta):=\nabla_{\bm{\xi}}\otimes\bm{x}$. This gives, at $\zeta=0$,

[TABLE]

By (4.1) we explicitly obtain

[TABLE]

so

[TABLE]

We can now introduce finite element spaces constructed from the basis previously discussed by defining

[TABLE]

We also need the set of interior edges defined by

[TABLE]

and the set of boundary edges on the Dirichlet part of the boundary

[TABLE]

To each interior edge $E$ we associate the conormals $\bm{\nu}^{\pm}_{E}$ given by the unique unit vector which is tangent to the surface element $K^{\pm}$, perpendicular to $E$ and points outwards with respect to $K^{\pm}$. Note that the conormals $\bm{\nu}_{E}^{\pm}$ may lie in different planes at junctions between different plates. The jump and average of multilinear forms for edges $E\in\mathcal{E}_{h}^{I}$ are defined by (3.12). For edges $E\in\mathcal{E}_{h}^{D}$ it is convenient to use the notation

[TABLE]

4.2 The Method

Our finite element method takes the form: find ${\bm{U}}\in\bm{V}_{h}:=[W_{h}]^{3}$ such that

[TABLE]

Here the bilinear form $A_{h}(\cdot,\cdot)$ is defined by

[TABLE]

with $\bm{v}=v_{n}\bm{n}+\bm{v}_{t}$ and

[TABLE]

where $(\cdot,\cdot)_{\omega}$ denotes the $L_{2}(\omega)$ scalar product, and

[TABLE]

Here $\beta=\beta_{0}(2\mu+2\lambda)$ where $\beta_{0}$ is an $O(1)$ constant, cf. HaLa11 , and we also recall that the factor $\tilde{t}^{2}$ is included in the definition (3.9) of the moment $M$. The right hand side is given by

[TABLE]

This is a c/dG method closely related to the one studied in HaLa11 , with the difference of being formulated in an arbitrary orientation in $\mathbb{R}^{3}$, including membrane deformations, and extended to structures of plates.

We note that:

• •

The continuity of displacement (3.22) is strongly enforced since $\bm{V}_{h}$ consists of continuous functions.

• •

The continuity of the rotation angle (3.23) is weakly enforced by the discontinuous Galerkin method.

• •

The force equilibrium conditions (3.24) and (3.26) are weakly enforced but does not give rise to any additional terms in the formulation since $\bm{V}_{h}$ consists of continuous functions.

• •

The moment equilibrium condition (3.25) is weakly enforced by the discontinuous Galerkin method.

More precisely, consider an edge $E\in\mathcal{E}^{I}_{h}$ shared by two elements $K^{+}$ and $K^{-}$. Multiplying the exact equation by a test function $\bm{v}\in\bm{V}_{h}$ and using Green’s formula element wise generates the following contribution at the edge $E$,

[TABLE]

where $\bm{F}^{\pm}=\bm{F}^{\pm}(\bm{u})$ and $M^{\pm}=M^{\pm}(\bm{u})$. For the first term we have using the continuity of $\bm{v}$ and (3.6),

[TABLE]

For the second term we note that the integrand may be written

[TABLE]

where we used the fact that $M^{\pm}$ is $2$-linear in $\bm{\nu}^{\pm}$, see (3.9), and $\bm{\nu}^{\pm}\cdot\nabla_{\Gamma}v^{\pm}_{n}$ is 1-linear in $\bm{\nu}$, and thus $M^{\pm}\bm{\nu}^{\pm}\cdot\nabla_{\Gamma}v^{\pm}_{n}$ is 3-linear in $\bm{\nu}^{\pm}$, together with the definition (3.12) of the jump to write the sum as a jump. Next using (3.16) we get

[TABLE]

since $[M]=0$ according to (3.25). Thus the second term takes the form

[TABLE]

where at last we symmetrized using the fact that the added term is zero by (3.23) and we included the dependency $M=M(\bm{u})$ for clarity. We finally note that we have the following identities

[TABLE]

and

[TABLE]

Remark 4.1

We note that the method for a plate structure has the same form as for a single plate since we use the proper definitions of jumps and averages encoded by the conormal.

Remark 4.2

We note that with this formulation, we have Galerkin orthogonality

[TABLE]

which enables us to prove an a priori error estimate of optimal order provided the solution is regular enough using the same techniques as in HaLa02 .

Remark 4.3

For shell modelling, the plate approach can still be used by viewing the shell as an assembly of facet elements. Then we have an elementwise planar approximation $\Gamma_{h}$ of $\Gamma$ and we use elementwise projections $\bm{P}_{h}=\bm{I}-\bm{n}_{h}\otimes\bm{n}_{h}$, where $\bm{n}_{h}$ is the elementwise constant approximation of $\bm{n}$. The differential operators are then defined on the discrete surface, e.g, $\nabla_{\Gamma_{h}}v:=\bm{P}_{h}\nabla v$, etc., and replacing the exact differential operators and exact surface by their discrete approximations in (4.12) we obtain a simple shell model.

5 Numerical Examples

We consider the surface of the box $[0,1]\times[0,1]\times[0,1]$, fixed to the floor and with one wall missing. The material data are: Poisson’s ratio $\nu=0.5$ and Young’s modulus $E=10^{9}$. The stabilization parameter was set to $\beta_{0}=10$. An ad hoc residual–based adaptive scheme was used to generate locally refined meshes. The load was given as

[TABLE]

at $x=0$, $\bm{f}=\bm{0}$ elsewhere. The point of the scaling with thickness is that after division by $t^{2}$ the membrane stiffness will scale with $t^{-2}$ so that the limit of $t\rightarrow 0$ corresponds to the inextensible plate solution. With increasing $t$ the membrane effect will become more and more visible. The numerical results using three different thicknesses, $t=10^{-k}$, $k=3,2,1$, are given in Figs. 1–3. Note the marked membrane deformations at $k=1$.

6 Concluding Remarks

In this paper we have introduced a c/dG method for arbitrarily oriented plate structures. Our method is expressed directly in the spatial coordinates, unlike traditional schemes that typically are based on coordinate transformations from planar elements. This leads to a remarkably simple and easy to implement discrete scheme. The c/dG approach also allows for avoiding the use of $C^{1}$–continuity, otherwise required by the plate model, by allowing for discontinuous rotations between elements, and the same function space can then be used to model both plate and membrane deformations. We also introduced the proper conormals, mean values, and jumps necessary for handling the discontinuities on the element borders.

Acknowledgements.

This research was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants Nos. 2011-4992, 2013-4708, and Swedish strategic research programme eSSENCE.