An ultraweak formulation of the Kirchhoff–Love plate bending model and DPG approximation Supported by CONICYT through FONDECYT projects 1150056, 11170050, The Magnus Ehrnrooth foundation, and by Oulun rakennustekniikan säätiö.

An ultraweak formulation of the Kirchhoff–Love plate bending model and DPG approximation thanks: Supported by CONICYT through FONDECYT projects 1150056, 11170050, The Magnus Ehrnrooth foundation, and by Oulun rakennustekniikan säätiö.

Thomas Führer    Norbert Heuer Facultad de Matemáticas, Pontificia Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, Santiago, Chile, email: {tofuhrer,nheuer}    Antti H. Niemi Structures and Construction Technology Research Unit, Faculty of Technology, University of Oulu, Erkki Koiso-Kanttilan katu 5, Linnanmaa, 90570 Oulu, Finland, email:

We develop and analyze an ultraweak variational formulation for a variant of the Kirchhoff–Love plate bending model. Based on this formulation, we introduce a discretization of the discontinuous Petrov–Galerkin type with optimal test functions (DPG). We prove well-posedness of the ultraweak formulation and quasi-optimal convergence of the DPG scheme.

The variational formulation and its analysis require tools that control traces and jumps in (standard Sobolev space of scalar functions) and (symmetric tensor functions with -components whose twice iterated divergence is in ), and their dualities. These tools are developed in two and three spatial dimensions. One specific result concerns localized traces in a dense subspace of . They are essential to construct basis functions for an approximation of .

To illustrate the theory we construct basis functions of the lowest order and perform numerical experiments for a smooth and a singular model solution. They confirm the expected convergence behavior of the DPG method both for uniform and adaptively refined meshes.

Key words: Kirchhoff–Love model, plate bending, biharmonic problem, fourth-order elliptic PDE, discontinuous Petrov–Galerkin method, optimal test functions

AMS Subject Classification: 74S05, 74K20, 35J35, 65N30, 35J67

1 Introduction

We develop an ultraweak variational formulation for a bending-moment variant of the Kirchhoff–Love plate model, and present a discontinuous Petrov–Galerkin method with optimal test functions (DPG method) that is based on this formulation. We prove well-posedness of the continuous formulation and quasi-optimal convergence of the discrete scheme. At the heart of the analysis is the space and its traces and jumps. This space consists of symmetric tensors with -components whose twice iterated divergence is in (the notation indicates the divergence operator that acts on the rows of tensors).

The Kirchhoff–Love model was introduced by Kirchhoff [32] in a form that is generally accepted today. Kirchhoff also applied the model to determine the free vibration frequencies and modes of circular plates. A historical account of the development of the model is incorporated in [33] where Love uses Kirchhoff’s approach to study vibrations of initially curved shells. Nowadays, the model is widely used in structural engineering, e.g., to dimension reinforced concrete slabs under static loads [28] and to control disturbing vibrations of wooden floors and other lightweight plane structures.

Perhaps the most well-known mathematical representation of the Kirchhoff–Love model for linearly elastic and isotropic material is given by the biharmonic equation

where is the deflection of the plate mid-surface , is the Laplace operator and and represent the external loading and bending rigidity of the plate, respectively.

It is evident that application of the model to complex geometries requires employment of numerical methods such as the finite element method. The literature on the numerical analysis of plate bending problems is vast due to the aforementioned practical relevance of the problems and respectable age of the structural models. It is not feasible to perform a thorough literature review here but two points that motivate our work can be made. First, conventional methods based on the variational principle of virtual displacements produce as direct output only the deflection values. These, albeit needed values, are not sufficient for structural design purposes where stresses and their resultants are of utmost importance. Second, verification of numerical accuracy of finite element algorithms is at the hearth of simulation governance, see [37]. This is a serious challenge in practical plate-bending problems where both the geometry and applied loading can be very irregular so that many of the contemporary developments in the finite element modeling of plate problems are devoted to a posteriori error estimation and adaptivity, see, e.g., [9, 2, 27].

We develop the theoretical framework for a DPG discretization to address the above challenges and, perhaps more significant, to set a theoretical basis to develop and analyze DPG schemes for other structural models like the singularly perturbed Reissner–Mindlin plate model and different shell models. Our analysis includes the case of singular problems on non-convex plates in contrast to many publications that assume convexity or smooth boundaries. In this context, we mention the mixed formulation from Amara et al. [1] who specifically use the space (without symmetry), thus allowing for singularities. Their numerical scheme is based on a decomposition of resulting in a mixed formulation that can be discretized by standard finite elements. In [25], Gallistl proposes a similar splitting approach for polyharmonic problems with corresponding finite element scheme.

The DPG framework has been founded by Demkowicz and Gopalakrishnan in [15]. It is very flexible and can be used with various variational formulations. A posteriori error estimation is also built-in, see [18]. DPG schemes have been applied previously to structural engineering problems in [34, 8] and to more general problems of elasticity in [3, 31]. The most closely related investigation to the present work is probably [8]. That investigation showed that an ultraweak variational formulation of the Reissner–Mindlin plate bending model is well posed and that the associated discretization is convergent. Rather accurate numerical results were observed despite the fact that the theoretically obtained stability constant is very weak and depends on the slenderness of the plate. In particular, the question of well-posedness of the ultraweak variational formulation of the asymptotic Kirchhoff–Love model was left open.

Essential motivation for the use of DPG schemes is their possible robustness for singularly perturbed problems. The intrinsic energy norm can bound in a robust way approximation errors in the sense that quasi-optimal error estimates (by the energy norm, which is accessible) are uniform with respect to perturbation parameters. To achieve this robustness in appropriate (selected) norms, it is paramount to have an appropriate variational formulation, and proving robustness is usually non-trivial. For an analysis of second-order elliptic problems with convection-dominated diffusion (“confusion”) and reaction-dominated diffusion (“refusion”) we refer to [19, 12, 5, 6] and [30], respectively. The DPG setting for refusion from [30] has been extended to transmission problems and the coupling with boundary elements [23], and to Signorini-type contact problems [24].

Whereas we do not consider a singularly perturbed problem in this paper, the development of a DPG scheme for the Kirchhoff–Love model is relevant in its own right as discussed before, and will be essential to deal with other models of plate problems. Since we expect our technical tools to be useful also for fourth-order problems in three dimensions, they are developed for both two and three space dimensions (they can be generalized to any space dimension). Discretizations of fourth-order problems usually avoid -bilinear forms to employ simpler than -conforming basis functions. In this respect, our choice of ultraweak variational formulation has the advantage that field variables are only in -spaces whereas appearing trace variables (traces of and ) are relatively straightforward to discretize.

Let us discuss the structure of our work. In the next section we introduce the model problem of a certain bending-moment formulation for the Kirchhoff–Love model. For simplicity we assume fully clamped plates but this is not essential as our formulation gives access to all kinds of boundary conditions. In that section, we also start developing a variational formulation. Since DPG schemes use product spaces111Often they are referred to as “broken” spaces. We prefer to call them product spaces since important Sobolev spaces, e.g., of negative order or order , cannot be localized but have to be defined as product spaces from the start, cf. [29]. with respect to subdivisions of into elements, trace operations in the underlying Sobolev spaces appear naturally. For fourth-order problems this is a non-trivial issue. Therefore, in order to define a well-posed variational formulation in product spaces we need to develop trace and jump operations, in our case in and . This is subject of Section 3, whose contents is discussed in more detail below. Eventually, in Section 4, we are able to define our ultraweak variational formulation and state its well-posedness (Theorem 11). We then briefly define the DPG scheme and state its quasi-optimal convergence (Theorem 12). Proofs of Theorems 11 and 12 are given in Section 5. We do not dwell much on the discussion of DPG schemes and their analysis. It is known that an analysis of the underlying adjoint problem gives access to the well-posedness of the variational formulation and quasi-optimal convergence of the DPG method (references have been given above). Though we do stress the fact that our analysis goes beyond standard techniques. Rather than splitting the adjoint problem into a homogeneous one in product spaces and an inhomogeneous one in global (“unbroken” or non-product) spaces (like, e.g., in [14, 19, 30]) or deducing stability of the adjoint problem in product spaces from the one of the global form [11], we consider the full adjoint problem as a whole. Section 5 starts with defining the adjoint problem. Its well-posedness is proved in §5.1. Key idea is to describe the primal unknown of the adjoint problem as the solution to a saddle point problem without Lagrange multiplier. Specifically, the primal unknown stays in the original product space and test functions are considered in the corresponding global space. Of course, this problem could be reformulated as a traditional saddle point problem. However, our technique is applicable to adjoint problems with data that require continuity,222Here we only note that such restrictions appear when considering first-order formulations of plate bending models. that is, leaving the setting of ultraweak formulations. In this sense, our new technique of analyzing the adjoint problem is fundamental. Extensions to other problems will be subject of future research.

Let us note that there is a recent abstract framework by Demkowicz et al. [17]. Under specific assumptions it yields the well-posedness of -ultraweak formulations in product spaces without explicitly analyzing trace spaces. In [26], Gopalakrishnan and Sepúlveda applied this setting to acoustic wave problems. In both references, an essential density assumption is only proved for simple geometries. Furthermore, trace variables are discretized via their domain counterparts whereas we only discretize the traces. It is also unknown whether the new framework gives robust control of variables in the case of singularly perturbed problems. In [21, 20], Ernesti and Wieners presented a simplified DPG analysis based on the framework from [17]. They use the density results for simple geometries from [17, 26]. Furthermore, the construction of their trace discretization is done without explicitly defining the domain parts, although they are needed for the stability and approximation analysis. In conclusion, in comparison with the current state of the framework from [17], our strategy has the advantages of giving access to singularly perturbed problems, being extendable to non- settings, avoiding domain contributions for trace discretizations, and not requiring density assumptions which can be hard to prove (though, see [1, Proposition 2.1] for the density of smooth tensor functions in defined by the graph norm without symmetry).

Now, to continue discussing the contents of our paper, having the analysis of the adjoint problem from §5.1 at hand, the proofs of Theorems 11 and 12 are straightforward. They are given in §5.2. Finally, in Section 6 we discuss the construction of discrete spaces for our DPG scheme and give some numerical examples. §6.1 is devoted to the construction of lowest-order basis functions. Whereas the field variables do not require any continuity across element interfaces, it is more technical to identify unknowns associated with trace variables. Specifically, the construction of basis functions for traces of requires to identify local continuity constraints. It turns out that traces of -functions cannot be split into natural components that allow for such a construction. This is analogous to where one uses a slightly more regular subspace of vector functions with normal (then localizable) traces in . In the literature, this subspace is usually denoted by . In the situation is worse since the definition of traces requires to integrate by parts twice. This generates two combined traces. We present lowest-order basis functions (for traces of ) that correspond to local unknowns associated with edges and nodes of triangular elements, plus jump constraints associated with interior nodes and neighboring elements. These constraints can be imposed by Lagrange multipliers. For sufficiently smooth solutions, our lowest order scheme converges with optimal order (Theorem 19). This result assumes the use of optimal test functions whereas, obviously, our numerical implementation uses approximated optimal test functions. We do not analyze the influence of this approximation here. In §6.2 we present numerical results for two examples, the case of a smooth solution and the case of a singular solution. Uniform mesh refinement yields optimal and sub-optimal convergence, respectively, whereas an adaptive variant restores optimal convergence for the singular example. It is worth mentioning that the singular example solution generates a tensor of whose divergence is less than -regular. This shows, in particular, that our analysis of traces and jumps in cannot be split into two steps/spaces (symmetric tensors in whose divergence are elements of ).

To conclude, the central focus of this paper is on the analysis of traces and jumps in , in Section 3. Despite of considering a plate model, this analysis is done in two and three space dimensions. It is relevant for other fourth-order problems in three dimensions. Section 3 is split into several subsections. In the first two, §§3.1 and 3.2, we define and analyze trace operators in and (denoted by and , with local versions and , respectively), and corresponding trace spaces and norms. In §3.3, we consider the product variant of and jumps of its elements. Specifically, we characterize the inclusion through (vanishing) duality with (Proposition 4). In §3.4 we revisit (a subspace of) the product space and study traces rather than jumps (of course, trace operators can be used to define and analyze jumps). We define a dense product subspace and prove that our previous “local” trace operators (they act on boundaries of elements) can be further localized when acting on this subspace (Proposition 6). This is of utmost importance for the numerical scheme since it implies density of our discrete spaces in , and thus convergence. §3.5 corresponds to §3.3, considering jumps of a product space rather than of , with continuity characterization by duality with the trace space (Proposition 8).

The final Subsection 3.6 provides a Poincaré inequality in the product space . Recall that traditional stability proofs of adjoint problems separate the analysis into a global non-homogeneous problem and a homogeneous one in product spaces and with jump data. The non-homogeneous problem usually gives control of a seminorm of the primal variable so that a Poincaré inequality is required to bound the norm. Furthermore, proving stability of homogeneous adjoint problems with jump data is usually done via a Helmholtz decomposition. For details see, e.g., [14, Lemmas 4.2, 4.3]. In our case, the global adjoint problem gives also only access to a seminorm of the primal variable, and still, the connection between jump data and the field variable is established by a Helmholtz decomposition. We combine both techniques and give a short proof of a Poincaré inequality in which uses a Helmholtz decomposition only implicitly.

Throughout the paper, means that with a generic constant that is independent of the underlying mesh (except for possible general restrictions like shape-regularity of elements). Similarly, we use the notation and .

2 Model problem

We start by recalling the Kirchhoff–Love model, cf. [38]. The static variables of the model are the shear force vector and the symmetric bending moment tensor . These stand for stress resultants representing internal forces and moments per unit length along the coordinate lines on the plate mid-surface . They are related to the external surface load and to each other by the laws of static equilibrium (force and moment balance) as

The operator denotes the divergence of vector functions, and is the divergence operator acting on rows of tensors. Denoting by the infinitesimal strain tensor, or symmetric gradient, we introduce the bending curvature , the Hessian of in our case. For linearly elastic isotropic material, the bending moments can be written in terms of as


is the bending rigidity of the plate defined in terms of the Young modulus and Poisson ratio of the material and the plate thickness . The values of these parameters are not very critical concerning the numerical solution of the problem. acts as scaling parameter and the influence of the Poisson ratio on the solution is mild. We select fixed and so that is positive definite.

Let us now assume that () is a bounded simply connected Lipschitz domain with boundary . (Of course, for the plate-bending problem, only is physically motivated.) For a given our model problem is


Here, is the exterior unit normal vector on . Later, will be used generically for normal vectors. Before starting to develop a variational formulation, we introduce a mesh that consists of general non-intersecting open Lipschitz elements. Only in §3.4 we will require that the mesh is conforming and consists of generalized (curved) polyhedra/polygons, and in the numerical section §6 we restrict ourselves to two space dimensions and conforming triangular meshes of shape-regular elements. To the mesh we associate the skeleton . For , scalar functions and symmetric tensors , let us define the norms

and induced spaces , as closures of and with respect to the corresponding norm. Here, and are the spaces of smooth functions and smooth symmetric tensors, respectively, on . (The logic for the notation with plain letter is that tensors are mapped to scalar functions by the operator . Similarly, below we introduce with bold as maps tensor functions to vector functions.) Throughout the paper, denotes the -norms for scalar, vector and tensor functions on the indicated set . When we drop the index and simply write instead of . The corresponding bilinear forms are and . The spaces and denote symmetric tensor functions on and , respectively.

Now, given a mesh , we define product spaces (tacitly identifying product spaces with their broken variants)

with canonical product norms and , respectively. We will also need the global spaces and which are the closures of and , respectively, with corresponding norms and , and similarly for .

Now, we test

Formally integrating by parts on every element and summing over the elements and summing the two equations, the testing results in


Here and in the following, a differential operator with index means that it is taken piecewise with respect to the elements . We will write equivalently, e.g., , and similarly for other differential operators taken in a piecewise form. Furthermore, we use the generic notation for the unit normal vector on and , pointing outside and , respectively. The notation , and later , indicate dualities on and , respectively, with -pivot space.

At this point it is not clear whether the appearing normal components in (2) on the boundaries of elements are well defined. Indeed, essential part of this paper is to study the relation between traces and jumps of the involved spaces , , and . This will be done in the next section, before returning to a variational formulation of (1) in Section 4.

3 Traces, jumps and a Poincaré inequality

In the following we introduce and analyze operators and norms that serve to give the terms , , , and , from (2) a meaning for and .

3.1 Traces of

We start by introducing linear operators for by


We note that this definition is consistent with the observation made by Amara at al. in [1, Theorem 2.2] (they consider the whole domain instead of an element ). The range of the operator is denoted by

These traces are supported on the boundary of the respective element since

cf. Proposition 4 below. It is therefore clear that, for given , the duality only depends on the traces of and on . Analogously, for any (smooth symmetric tensors with support in ) since

Remark 1.

Let us define as the closure of with respect to the norm . It is clear that is not a subspace of , nor is a subspace of (see the second example in §6.2). This is precisely the reason we have to consider the trace operator in the form (3). When restricting this operator as

it reduces to standard trace operations. In this case the two dualities are defined independently in the standard way,

Now, in the three-dimensional case , defining the tangential trace for and the surface gradient , we formally write


Correspondingly, in two dimensions (), we introduce the unit tangential vector along in mathematically positive orientation, and use the notation for with corresponding tangential derivative . Then, (4) applies as well. We also need the surface divergence defined by for sufficiently smooth vector functions with . For precise definitions and appropriate spaces we refer to [7]. With these definitions it is clear that we can define separate traces




that coincide with the corresponding trace terms for sufficiently smooth functions , cf. the operators and in [1, page 1635]. On the one hand, these traces are relevant to identify basis functions for the approximation of traces of -functions. On the other hand, specifying one of these traces, the other is well defined as a functional acting on traces of -functions (without the trace conditions for in (5) and (6)). Applying this on the boundary of , it is possible to specify any physically meaningful boundary condition based on the terms , , , and on . Note that refers to the operator that is dual to the (negative) global surface gradient , and integrating by parts piecewise on subsets of generates a piecewise surface divergence plus jumps at the interfaces, cf. [35]. Indeed, these jumps will be essential for the approximation analysis, based on Proposition 6 below.

The collective trace operator is defined by

with duality


and range

(Here, and in the following, considering dualities on the whole of , possibly involved traces onto are always taken from without further notice.) These global and local traces are measured in the minimum energy extension norms,

Alternative norms are defined by duality as follows,


Here, for given , the duality with is defined by

and, for and ,


This is consistent with definitions (3) and (7).

Lemma 2.

It holds the identity

so that

has unit norm and is closed.


The estimate is immediate by bounding

Now let and be given. We define as the solution to the problem


Note that the right-hand side functional implies a natural boundary condition for . Furthermore, since for , satisfies


first in the distributional sense and, by the regularity of , also in . Using the function we continue to define as the solution to


Again, the right-hand side functional induces a natural boundary condition for , and it holds


We show that . Indeed, defining , we find with (11) that so that by definition of and definition (3) of ,

for any . This shows that solves (12) and by uniqueness, . Using this relation and , it follows by (10) that

In other words, . This relation together with selecting in (10) and in (12), shows that


Noting that

by (13), relation (14) finishes the proof of the norm identity. The space is closed as the image of a bounded below operator. ∎

3.2 Traces of

Let us study traces of in a similar way as in the previous section.

We define linear operators for by


and observe that (cf. (3))


The ranges are denoted by

It is immediate that if and only if . The collective trace operator is defined by

with duality


and range

These trace spaces are provided with canonical trace norms,

Alternative norms are defined by duality as follows,

Here, for given , the duality with is defined by

and, for and ,


This is consistent with definitions (15) and (17).

Lemma 3.

It holds the identity

so that

has unit norm and is closed.


The proof is very similar to that of Lemma 2.

The estimate follows by bounding

To show the other inequality, let and be given. We define as the solution of


It satisfies


We continue to define by


It satisfies


and we conclude that as follows. Defining , (20) shows that so that by definition of and (cf. (15)),

for any . Hence solves (21), that is, . Using this relation and , (19) shows that

so that . Then selecting in (19) and in (21), we obtain



by (22), relation (23) finishes the proof of the norm identity. The space is closed as the image of a bounded below operator. ∎

3.3 Jumps of

Proposition 4.

(i) For it holds

(ii) The identity

holds true.


The proof of (i) follows the standard procedure, cf. [11, Proof of Theorem 2.3]. For and ,

showing the direction “”. Now, for given with for any we have in the distributional sense

Therefore, , that is, .

Next we show (ii). The inequality holds by definition of the norms. To show the other inequality let be given. By definition of there is such that . Furthermore, for , let be such that and

Then, defined by () satisfies (cf. (16))