A simple approach for finite element simulation of reinforced plates

We present a new approach for adding Bernoulli beam reinforcements to Kirchhoff plates. The plate is discretised using a continuous/discontinuous finite element method based on standard continuous piecewise polynomial finite element spaces. The beams are discretised by the CutFEM technique of letting the basis functions of the plate represent also the beams which are allowed to pass through the plate elements. This allows for a fast and easy way of assessing where the plate should be supported, for instance, in an optimization loop.


Introduction
Reinforcements of solids using lower-dimensional structures such as beams can be simulated in a finite element context by coupling the variables of the beam to the variables of the solid, either along element edges as in McCune, Armstrong, and Robinson [13] or by interpolation on element edges as in Sadek and Shahrour [15]. In the latter case, the beam geometry can be modelled independently of the bulk mesh which is crucial; however, the finite element approximation of the lower-dimensional object is otherwise independent and uncoupled to the solid, and the rotation degrees of freedom of beams are hard to match to the solid (if they are to influence the solution in the solid).
In [11] we proposed to use the same finite element space for the beam as for the higher dimensional structure; more precisely, the trial and test space for the beam is obtained by taking the restriction or trace to the beam. Here we further develop this approach to allow for coupling between plates and beams, more precisely the Kirchhoff-Love plate model and the Bernoulli beam model. These models involve fourth order partial differential equations. We discretize these models using the so called continuous/discontinuous Galerkin, c/dG, method which relaxes the required C 1 continuity of the shape functions for the beam and plate by use of a discontinuous Galerkin approach with C 0 -continuity. We emphasise that the concept is quite general, as illustrated in our previous work on embedding in elastic  solids, of membranes [4] and of embedded trusses and beams [11]. A similar approach was recently suggested for modelling embedded trusses by Lé, Legrain, and Moës [12].
2. Modeling of Reinforced Plates 2.1. The Basic Approach. In this Section we develop a simple model of a set of beam elements in a plate. The main approach is as follows: • Given a continuous finite element space, based on at least second order polynomials for the plate, we define the finite element space for the one-dimensional structure as the restriction of the plate finite element space to the structure which is geometrically modeled by an embedded curve or line. • To formulate a finite element method on the restricted or trace finite element space we employ continuous/discontinuous Galerkin approximations of the Euler-Bernoulli beam model. The beams are then modeled using the CutFEM paradigm and the stiffness of the embedded beams is in the most basic version, which we consider here, simply added to the plate stiffness. To ensure coercivity of the cut beam model we in general need to add a certain stabilization term which provides control of the discrete functions variation in the vicinity of the beam. However, for beams embedded in a plate, the plate stabilizes the beam discretizations, and we shall show that if the plate is stiff enough compared to the beam the usual additional stabilization [1] is superfluous. The plate problem may also be viewed as an interface problem in order to more accurately approximate the plate in the vicinity of the beam structure; this approach is however significantly more demanding from an implementation point of view and we leave it for future work.
The work presented here is an extension of earlier work [4] where membrane structures were considered, in which case a linear approximation in the bulk suffices.
2.2. The Kirchhoff-Love Plate Model. In the Kirchhoff-Love plate model, posed on a polygonal domain Ω ⊂ R 2 with boundary ∂Ω and exterior unit normal n, we seek an out-of-plane (scalar) displacement u to which we associate the strain (curvature) tensor and the plate stress (moment) tensor 12(1 + ν Ω ) with E Ω the Young's modulus, ν Ω the Poisson's ratio, and t Ω denotes the plate thickness. Since 0 ≤ ν Ω ≤ 0.5 the constants are uniformly bounded.

2.3.
The Euler-Bernoulli Beam Model. Consider a straight thin beam with centerline Σ ⊂ Ω and a rectangular cross-section with width b Σ and thickness t Σ , see Figure 2. The modeling of the beam is performed using tangential differential calculus and we follow the exposition in [10,11], which also covers curved beams. Using this approach the beam equation is expressed in the same coordinate system as the plate, which is convenient in the construction of the cut finite element method for reinforced plates.
Let t be the tangent vector to the line Σ, embedded in R 2 . We let p : R 2 → Σ be the closest point mapping, i.e. p(x) = y where y ∈ Σ minimizes the Euclidean norm |x − y| R 3 . We define ζ as the signed distance function ζ(x) := ±|x − p(x)|, positive on one side of Σ and negative on the other.
Let P Σ = t ⊗ t be the projection onto the one dimensional tangent space of Σ and define the tangential derivatives Then we have the identity Based on the assumption that planar cross sections orthogonal to the midline remain plane after deformation we assume that the displacement takes the form where θ : Σ → R is an angle representing an infinitesimal rotation, assumed constant in the normal plane. In Euler-Bernoulli beam theory the beam cross-section is assumed plane and orthogonal to the beam midline after deformation and no shear deformations occur, which means that we have This definition for θ in combination with (13) constitutes the Euler-Bernoulli kinematic assumption We assume the usual Hooke's law for one dimensional structural members where E Σ is the Young modulus and the tangential strain tensor is given by where in the last equality we used the identity Next note that the strain energy density can be written and the total energy of the beam structure is obtained by integrating over the beam volume where the integral over the cross section is accounted for by the cross-section area and its second moment We are thus led to introducing the beam stress tensor and thus we have the beam Hooke law Taking variations we obtain the weak statement, assuming zero displacements and rotations at the end points of Σ, we thus seek where the forms are defined by Here we recognize the right hand side as the traditional bilinear form associated with the Euler-Bernoulli beam.
Remark 2. We note that in the alternative reinforcement geometry, right in Figure 2, we have We may also consider more complicated cross sections and compute the proper parameters.
2.4. The Reinforced Plate Model. Let S = {S } be a set of beams arbitrarily oriented in Ω. Using superposition we obtain the problem: find u ∈ V such that and the forms are defined by Remark 3. Note that for the alternative plate reinforcement geometry, right in Figure 2, there is no geometric error in our method if we use the parameters (30). In the standard reinforcement geometry, left in Figure 2, there is a however a geometric error proportional to b Σ in the plate bilinear form, which arises in the superposition since the intersection between the beam and the plate is nonempty. We will later see that b Σ typically is smaller (in practice significantly smaller) than the mesh size since we are using thin beam and plate theory, see (35), and thus the geometric error is small. • We consider a subdivision T h = {T } of Ω into a geometrically conforming finite element mesh, with mesh parameter h ∈ (0, h 0 ]. We assume that the elements are shape regular, i.e., the quotient of the diameter of the smallest circumscribed sphere and the largest inscribed sphere is uniformly bounded. We denote by h T the diameter of element T and by h = max T ∈T h h T the global mesh size parameter. • Since we are using thin plate and beam theory we assume that there is a constant C mesh such that • We shall use continuous, piecewise polynomial approximations, for both the membrane and plate problem. Let where P k (T ) is the space of polynomials of degree less or equal to k defined on T . For simplicity, we write V Ω,h = V Ω,h,k . • To define our method we introduce the set of faces (edges) F in the mesh, F h = {F}, and we split F h into two disjoint subsets where F h,I is the set of faces in the interior of Ω and F h,B is the set of faces on the boundary. • Further, with each face F we associate a fixed unit normal n F such that for faces on the boundary n F is the exterior unit normal. We denote the jump of a function v at a face F by • The intersection points between Σ and element faces in F h (Σ) is denoted and we assume that this is a discrete set of points (thus excluding the case where any F ∈ F h coincides with a part of Σ).

3.2.
The c/dG Method for the Plate. We approximate the solution to the plate problem using the continuous/discontinuous Galerkin (c/dG) method: Find u h ∈ V Ω,h , with k ≥ 2, such that The bilinear form a Ω,h (·, ·) is defined by Here β Ω is a positive parameter of the form where β Ω,0 is a constant depending on the polynomial order k, see [9] for details, and h F is defined on each face F by with |T | the area of T and |F| the length of F.
Remark 5. Other boundary conditions for plates, for instance simply supported and free, can easily be included in the c/dG finite element method, see [7] for details.

Remark 6.
For v ∈ V Ω,h we have [∇v] = [n F · ∇v]n F since v is continuous across a face and v = 0 on ∂Ω. Therefore for all v, w ∈ V Ω,h , and we note that (n F · σ(∇u) · n F )| F is the bending moment at the edge F.

3.3.
The Cut c/dG Method for a Beam. We propose the following cut c/dG method.
the penalty parameter takes the form with β Σ,0 a parameter that only depends on the polynomial order, and s h , with positive parameters γ Σ,i , is a stabilization term which is added to ensure coercivity and stability of the stiffness matrix, cf. [1].
Remark 7. Using the identities t ⊗ t we note that a Σ,h can alternatively be written in the form which is the form in [5].
Remark 8. The terms on the discrete set P h (Σ) are associated with the work of the end moments on the end rotation which occur due to the lack of C 1 (Ω) continuity of the approximation, as in the plate model. See Remark 6.
Remark 9. We note that due to the stabilization this method works for a single beam, i.e. without being embedded in a plate. The basic principle is the same as for the trace finite element method proposed in [14] and the stabilized version proposed in [2]. When the beam is embedded in a plate, which is the case in this work, the need for the stabilization term is mitigated, and if the plate is sufficiently stiff we may omit the stabilization term, see Section 3.5 for further details.

The c/dG Method for the Reinforced Plate Model.
Recall that S = {Σ} is a set of beams arbitrarily oriented in Ω. Using superposition we obtain the problem: find u h ∈ V Ω,h such that where the forms are defined by 3.5. Coercivity for Reinforced Plates. In this section we study the coercivity of the c/dG method for the reinforced plate. We shall use the stability provided by the plate to prove stability of the reinforced plate, without the need of the stabilizing terms (γ Σ,1 = γ Σ,2 = 0). This is only possible as long as the mesh size h is larger than or equal to the beam with b Σ . When this condition is not satisfied, stability uniform in h is achieved only when the stabilizing terms are included (γ Σ,1 , γ Σ,2 > 0), using similar ideas as in [2,3]. Coercivity of the Plate. We first recall that the c/dG method for the plate is coercive. Introducing the energy norm there is a constant m P > 0 such that for β Ω large enough. Coercivity of the Reinforced Plate. Next turning to the reinforced plate we introduce the energy norm associated with the beam Then there is a constant m such that for β Ω and β Σ large enough. Verification of (57). Using the following two inequalities, which we verify below, for some constant C 1 > 0, and for β Σ large enough, we have where m = min(m P /3, C 1 m P /3, 1). Verification of (58). We note that, for x ∈ Σ ∩ T , T ∈ T h , we have the inverse inequality and thus we have the estimate We note, using the definitions (4) and (24) of C P and C B , that where we used the condition that the beam width is smaller than the mesh size (35) and thus the right hand side is a positive constant independent of the mesh size and so is C 1 . Verification of (59). First we have the equality To estimate we employ the inverse inequality (65) as follows where we used the inequality ab ≤ (δa 2 + δ −1 b 2 )/2 for δ > 0. We then obtain (59) as follows Here we choose: δ small enough to guarantee that where as above, see (69), C 1 > 0 independent of the mesh parameter h, and β Σ such that

Numerical Examples
In this Section, we give some elementary examples of what can be achieved with the presented technique. In all numerical examples we use polynomial order 2, f Σ = 0, and corresponding to the solution u = x 2 (1 − x) 2 y 2 (1 − y) 2 for a clamped plate unsupported by beams.
In order to handle more general boundary conditions we in particular need to be able to impose end displacements on the beam in the case of a free plate (we note that strongly imposed boundary conditions on the plate are also enforced on the beam). Zero displacement of the beam endpoints x E are imposed by adding penalty terms (50), whereβ Σ,0 is a penalty parameter. These terms suffice for optimal order convergence (of the beam approximation) in the case of second degree polynomial approximations since the shear forces required for energy consistency are third derivatives of displacements, and thus equal zero.

4.1.
Simply supported plate using beams with different supports. We consider a simply supported plate on the domain Ω = (0, 1) × (0, 1) with Young's modulus E Ω = 100, Poisson's ratio ν Ω = 1/2, and thickness t Ω = 0.1. The plate is supported by two beams oriented as in Fig. 4, one at x = 0.499 and one at at y = 0.499 (to avoid intersection with the mesh lines). The computational mesh is shown in Fig. 5 and in Fig.6 whe show a close-up of the intersection between the beams and the mesh.
For this problem we test two different supports for the beams: simply supported and fixed, and two different stiffnesses for the beams: E Σ = 100E Ω and E Σ = 1000E Ω . The thickness and width of the beam are equal and the same as the thickness of the plate. In Fig.  7 we show the results using E Σ = 100E Ω , with simply supported and fixed supports; in Fig.  8 we give the corresponding isolines, and in Fig. 9 we show the results using E Σ = 1000E Ω , with simply supported and fixed supports; in Fig. 10 we give the corresponding isolines.

4.2.
Plate only supported by beams. Next, we consider a plate with free boundaries, supported only by beams. All data for the plate are the same as in the previous example. The plate is supported by four beams positioned at 1/3 and 2/3 from each boundary as indicated in Fig. 11. The beams have the same dimension as previously, with Young's modulus E Σ = 100E Ω . The computational mesh is unstructured and shown in Fig. 12.
We first consider the case when the beams are clamped at x = 1 and free elsewhere. In Fig. 13 we see the corresponding deformation in elevation and isoline plot. Next we consider the case when all beams are clamped, Fig. 14, and simply supported, Fig.15. Note the the slight increase in central displacement for the latter.

Conclusions
We have formulated a continuous/discontinuous Galerkin method for beam reinforced thin plates. The method has the advantage that we can discretize both the beam and plate problem with the same standard finite element spaces of continuous piecewise polynomials defined on triangles (or quadrilaterals).