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

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.


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 [13] 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 [7], 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. [4], and further developed for plate models by Hansbo et al. [5,6,8,9] and by Wells and Dung [14].See also Larsson and Larson [12] 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 [10] 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 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 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 [1] 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 [11], 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 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.

Tangential Differential Calculus
Let Γ be a piecewise planar two-dimensional surface imbedded in R 3 , with piecewise constant unit normal n and boundary ∂Γ, split into a Neumann part ∂Γ N where forces and moments are known, and a Dirichlet part ∂Γ D where rotations and displacements are known.For ease of presentation we shall assume that ∂Γ N = ∅ and that we have zero displacements and rotations on the boundary.The case of ∂Γ N ∅ 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 Γ by ζ(x), for x ∈ R 3 , fulfilling ∇ζ = n, we can define the domain occupied by the shell by where t is the thickness of the shell, which for simplicity will be assumed constant.The closest point projection p : Ω t → Γ is given by the Jacobian matrix of which is where I is the identity and ⊗ denotes exterior product.The corresponding linear projector P Γ = P Γ (x), onto the tangent plane of Γ at x ∈ Γ, is given by and we can then define the surface gradient ∇ Γ as The surface gradient thus has three components, which we shall denote by For a vector valued function v(x), we define the tangential Jacobian matrix as and the surface divergence ∇ Γ • v := tr v ⊗ ∇ Γ .

Displacement and Strain
Upon loading, each point x ∈ Ω t , in the plate undergoes a displacement where u 0 and w are vector fields defined on Γ, u 0 arbitrary and w a tangential vector, w • n = 0 on Γ, or w = P Γ θ with θ arbitrary.Thus, neglecting in-plane extensions for the moment, we can write We introduce the strain tensor ε as and define the symmetric part of the tangential Jacobian as The in-plane strain tensor ε Γ is implemented using the following identity If we write where is the curvature tensor, cf.[2,3].For planar Γ, n is constant, and this simplifies to The total in-plane strain tensor is thus given by In [2,3] it is also shown that the mid-plane rotation in the absence of shear deformation is given by 2e Γ (u n n) • n, and for shear deformable inextensible shells we thus have the shear deformation vector It is is also easy to verify that In the tangential setting, the Kirchhoff assumption of zero shear deformations can therefore be written Furthermore, we find that and for inextensible plates we get and in this case we thus only obtain contributions to the strain energy from the displacement field

Variational Formulations
We shall assume isotropic stress-strain relations, where σ is the stress tensor, and plane stress conditions, for which the Lamé parameters λ and µ are related to Young's modulus E and Poisson's ratio ν via For the in-plane stress tensor we find, by projecting (2.29) from left and right, The potential energy of the plate is postulated as where σ : ε = i j σ i j ε i j for second order Cartesian tensors σ and ε.Integrating in ζ, we obtain Under the assumption of clamped boundary conditions, the corresponding variational problem is to find for all v ∈ H 2 0 (Γ).Introducing also membrane deformations, the total potential energy E tot of the plate must take into account both the bending energy E P and the membrane energy E M , so that E tot = E P + E M , where Since we wish to use a 3D Cartesian vector field we redefine u := u 0 and u n := n • u, make use of (2.17), and introduce the function space We are then led to the variational problem of finding u ∈ V such that we may write (2.37) in the more compact form t2 For implementation purposes we note that for n constant (2.41)

Strong Form
The corresponding strong form of the problem is to find u = and 3 Plate Structures

Forces and Moments
Consider first a subdomain polygonal subdomain ω ⊂ Γ of the plate Γ with boundary ∂ω consisting of line segments γ i .
Using Greens formula on ω we obtain where we used the identity v n = v • n and moved the normal to the first slot in the bilinear form.Letting τ be a unit tangent vector to ∂ω, we may split the last term on the right hand side of (3.2) in normal and tangent contributions as follows 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 γ i , with unit tangent and normal where ∂γ i consists of the two end points of the line segment γ i .We introduce the following notation for the normal and tangent components of the force and the moment at each of the line segments γ on ∂ω.Furthermore, we introduce the corner, or Kirchhoff, forces at a corner x associated with a line segment γ i , which has x as one of its endpoints and τ i is the unit tangent vector to γ i directed into x.We then have the identity where X(∂ω) is the set of corners on the polygonal boundary ∂ω and I(x) is an enumeration of the two linesegments that has x as one of its endpoints.

Jumps and Averages
Consider a line segment γ shared by two plates Γ + and Γ − .We note that the force F ± is an R 3 valued 1-form in ν ± and the moment M ± is an R valued 2-form in ν ± .More generally let w ± = w ± (ν ± , . . .ν ± ) be an R n valued d-linear form in ν ± .Then we define the jump and average at γ by Note that when both plates Γ + and Γ − reside in the same plane ν − = −ν + and we recover, using linearity and the simplified notation w ± (ν ± , . . ., ν ± ) = w(ν ± ), the standard jump and similarly for the average.Finally, let w ± i be an R n valued d i -linear form in ν ± , then we note that ( -linear form in ν ± and we have the identity where for n = 1 the scalar product is just usual multiplication of scalars.We may verify (3.16) by (3.20) (3.21)

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 γ where two plates Γ + and Γ − intersect we have the interface conditions corresponding to continuity of displacements, continuity of the rotation angle, equilibrium of forces, and equilibrium of moments.
Furthermore, at each corner x, not residing on the boundary of the structure, we require equilibrium of the Kirchhoff forces 0 = i∈I(x) where I(x) is an enumeration of the line segments that meet in the corner x and F ± x,i is the Kirchhoff force emanating from plate Γ ± i , 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.

The Mesh and Finite Element Space
Let K ⊂ R 2 be a reference triangle and let P 2 ( K) be the space of polynomials of order less or equal to 2 defined on K. Let Γ be triangulated with quasi uniform triangulation K h and mesh parameter h ∈ (0, h 0 ] such that each triangle K = F K ( K) is planar (a subparametric formulation).We let E h denote the set of edges in the triangulation.
We here extend the discontinuous Galerkin method of Dedner et al. [1] for the Laplace-Beltrami operator to the case of the plate.We recall that Γ is piecewise planar and thus n is a piecewise constant exterior unit normal to Γ.
For the parametrization of Γ we wish to define a map from a reference triangle K defined in a local coordinate system (ξ, η) to any given triangle K on Γ.Thus the coordinates of the discrete surface are functions of the reference coordinates inside each element, x Γ = x Γ (ξ, η).For any given parametrization, we can extend it to Ω t by defining where −t/2 ≤ ζ ≤ t/2 and n is the normal to Γ.We consider in particular a finite element parametrization of Γ as where x i are the physical location of the (geometry representing) nodes on the initial midsurface and ψ i (ξ, η) are affine finite element shape functions on the reference element.(This parametrization is of course exact in the case of a piecewise planar Γ.)For the approximation of the displacement, we use a constant extension, where u i are the nodal displacements, and ϕ i are piecewise quadratic shape functions.We employ the usual finite element approximation of the physical derivatives of the chosen basis {ϕ i } on the surface, at (ξ, η), in matrix representation, as where By (4.1) we explicitly obtain so We can now introduce finite element spaces constructed from the basis previously discussed by defining We also need the set of interior edges defined by and the set of boundary edges on the Dirichlet part of the boundary To each interior edge E we associate the conormals ν ± E given by the unique unit vector which is tangent to the surface element K ± , perpendicular to E and points outwards with respect to K ± .Note that the conormals ν ± E may lie in different planes at junctions between different plates.The jump and average of multilinear forms for edges E ∈ E I h are defined by (3.12).For edges E ∈ E D h it is convenient to use the notation

The Method
Our finite element method takes the form: find Here the bilinear form A h (•, •) is defined by where (•, •) ω denotes the L 2 (ω) scalar product, and Here β = β 0 (2µ + 2λ) where β 0 is an O(1) constant, cf.[9], and we also recall that the factor t2 is included in the definition (3.9) of the moment M. The right hand side is given by This is a c/dG method closely related to the one studied in [9], with the difference of being formulated in an arbitrary orientation in R 3 , including membrane deformations, and extended to structures of plates.
We note that: -The continuity of displacement (3.22) is strongly enforced since 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 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 ∈ E I h shared by two elements K + and K − .Multiplying the exact equation by a test function v ∈ V h and using Green's formula element wise generates the following contribution at the edge E, where F ± = F ± (u) and M ± = M ± (u).For the first term we have using the continuity of v and (3.6), For the second term we note that the integrand may be written where we used the fact that M ± is 2-linear in ν ± , see (3.9), and ν ± • ∇ Γ v ± n is 1-linear in ν, and thus M ± ν ± • ∇ Γ v ± n is 3linear in ν ± , together with the definition (3.12) of the jump to write the sum as a jump.Next using (3.16) we get since [M] = 0 according to (3.25).Thus the second term takes the form where at last we symmetrized using the fact that the added term is zero by (3.23) and we included the dependency M = M(u) for clarity.We finally note that we have the following identities Remark 41 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 42 We note that with this formulation, we have Galerkin orthogonality 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 [7].
Remark 43 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 Γ h of Γ and we use elementwise projections P h = I−n h ⊗n h , where n h is the elementwise constant approximation of n.The differential operators are then defined on the discrete surface, e.g, ∇ Γ h v := P h ∇v, etc., and replacing the exact differential operators and exact surface by their discrete approximations in (4.12) we obtain a simple shell model.

Numerical Examples
We consider the surface of the box [0, 1]×[0, 1]×[0, 1], fixed to the floor and with one wall missing.The material data are: Poisson's ratio ν = 0.5 and Young's modulus E = 10 9 .The stabilization parameter was set to β 0 = 10.An ad hoc residual-based adaptive scheme was used to generate locally refined meshes.The load was given as at x = 0, f = 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 → 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.

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.