Morley Finite Element Method for the von K\'{a}rm\'{a}n Obstacle Problem

This paper focusses on the von K\'{a}rm\'{a}n equations for the moderately large deformation of a very thin plate with the convex obstacle constraint leading to a coupled system of semilinear fourth-order obstacle problem and motivates its nonconforming Morley finite element approximation. The first part establishes the well-posedness of the von K\'{a}rm\'{a}n obstacle problem and also discusses the uniqueness of the solution under an a priori and an a posteriori smallness condition on the data. The second part of the article discusses the regularity result of Frehse from 1971 and combines it with the regularity of the solution on a polygonal domain. The third part of the article shows an a priori error estimate for optimal convergence rates for the Morley finite element approximation to the von K\'{a}rm\'{a}n obstacle problem for small data. The article concludes with numerical results that illustrates the requirement of smallness assumption on the data for optimal convergence rate.


Introduction
Short history of related work. The von Kármán equations [17] model the bending of very thin elastic plates through a system of fourth-order semi-linear elliptic equations; cf. [2,17,23] and references therein for the existence of solutions, regularity, and bifurcation phenomena. The papers [5,9,12,13,25,27,31,32] study the approximation and error bounds for regular solutions to von Kármán equations using conforming, mixed, hybrid, Morley, C 0 interior penalty and discontinuous Galerkin finite element methods (FEMs).
The obstacle problem is a prototypical example for a variational inequality and arises in contact mechanics, option pricing, and fluid flow problems. The location of the free boundary is not known a priori and forms a part of the solution procedure. For the theoretical and numerical aspects of variational inequalities, see [20,22]. A unified convergence analysis for the fourth-order linear two-sided obstacle problem of clamped Kirchhoff plates in [6][7][8] studies C 1 FEMs, C 0 interior penalty methods, and classical nonconforming FEMs on convex domains and, analyse the C 0 interior penalty and the Morley FEM on polygonal domains.
The obstacle problem for von Kármán equations with a nonlinearity together with a free boundary offers additional difficulties. The obstacle problem in [26,28,33] concerns a different plate model with continuation, spectral, and complementarity methods, while the papers [29,30] study conforming penalty FEM.

INTRODUCTION
The present paper is the first on the fourth-order semilinear obstacle problem of a (very thin) von Kármán plate. The article derives existence, uniqueness (under smallness assumption on data) and regularity results of the von Kármán obstacle problem. Nonconforming FEMs appear to be more attractive than the classical C 1 conforming FEMs, so this article suggests the Morley FEM to approximate the von Kármán obstacle problem and derives an optimal order a priori error estimate with the best approximation plus a linear perturbation.

Results and overview.
A smallness assumption on the data is derived in Section 2 to show that (1.2) is well-posed. The regularity results of Section 3 establish that any solution (u, v) to (1.2) satisfies u, v ∈ H 2 0 (Ω) ∩ H 2+α (Ω) ∩ C 2 (Ω) for the index 1/2 < α ≤ 1 with α = min{α ′ , 1} and the index α ′ of elliptic regularity [4] of the biharmonic operator in a polygonal domain Ω. Section 4 introduces the Morley FEM and discusses the well-posedness of the discrete problem with an a priori and an a posteriori smallness condition on the data for global uniqueness. Section 5 derives a priori energy norm estimates of optimal order α for the Morley FEM under the smallness assumption on the data that guarantees global uniqueness of the minimizer on the continuous level. The article concludes with numerical results that illustrates the requirement of smallness assumption on the data for optimal convergence rate.

Well-posedness
This section establishes the well-posedness of the problem (1.2). The existence of a solution to (1.2) follows with the direct method in the calculus of variations. The subsequent bound applies often in this paper and is based on Sobolev embedding. Let C S denote the Sobolev constant in the Sobolev embedding H 2 0 (Ω) ֒→ C(Ω) and let C F denote the Friedrichs constant with This means G(ξ) is the Riesz representation of the linear bounded functional b(ξ, ξ, •) in the Hilbert space (H 2 0 (Ω), a(•, •)). Consider the minimizer u of the functional j(ξ) for ξ ∈ K and 3) The equivalence of (1.2) with (2.3), for K = H 2 0 (Ω), is established in [17,Theorem 5.8.3]. Analogous arguments also establish the equivalence, for any non-empty, closed, and convex subset K of H 2 0 (Ω), so the proof is omitted. This implies that, to prove the existence of a solution to (1.2), it is sufficient to prove the existence of a minimizer to (2.3).
Proof. Given ξ ∈ K, the definition of j(•) in (2.3) and the Cauchy-Schwarz inequality lead to This implies the lower bound Consequently, there exists a sequence (u n ) n∈N in K such that The Cauchy-Schwarz and the Young inequalities lead to Consequently, |||u n ||| 2 + 2|||G(u n )||| 2 ≤ 4 j(u n ) + 4 f 2 −2 . Since j(u n ) is convergent, the sequences (u n ) n∈N and (G(u n )) n∈N are bounded in H 2 0 (Ω). Hence, there exist u, w ∈ H 2 0 (Ω) and a weakly convergent subsequence (u n k ) k ∈N such that u n k ⇀ u and G(u n k ) ⇀ w weakly in H 2 0 (Ω) as k → ∞.
The non-empty closed convex set K of H 2 0 (Ω) is sequentially weakly closed and so u ∈ K. Since u n k converges weakly to u in H 2 0 (Ω), this implies ∫ The compact embedding of H 2 0 (Ω) in L 2 (Ω) implies u n k → u in L 2 (Ω). Further for a given ϕ ∈ D(Ω), the definition of G(•) in (2.2), the symmetry of b(•, •, •) with respect to second and third arguments, and the weak convergence of u n k ⇀ u in H 2 0 (Ω) lead to Since ϕ is arbitrary in the dense set D(Ω) of H 2 0 (Ω), this means G(u n k ) ⇀ G(u) weakly in H 2 0 (Ω) as k → ∞. The sequentially weak lower semi-continuity of the norm ||| • ||| shows j(u) ≤ lim inf k j(u n k ). This and lim k→∞ j(u n k ) = β ≤ lim inf k j(u n k ) prove that u minimizes j in K. By the definition of G(•) in (2.2), (u, G(u)) solves (1.2). This concludes the proof. Theorem 2.3 establishes an a priori bound and the uniqueness of the solution to (1.2). Recall the Sobolev (resp. Friedrichs) constant C S (resp. C F ) from (2.1).

Theorem 2.3 (a priori bound and uniqueness)
.
Remark 2.1 (a priori and a posteriori criteria for uniqueness). The a priori smallness assumption on data C S M( f , χ) < √ 2 − 1 implies C S N( f , u) < √ 2 − 1 and so global uniqueness of the solution to (1.2). The first condition is a priori, but given the constant √ 2 − 1, f , and χ, M( f , χ) is hard to quantify. The second condition C S N( f , u) < √ 2 − 1 is a posteriori in the sense that u can be replaced by some ϕ ∈ K. Once an approximation u M of u is known, some ϕ ∈ K can be postprocessed by u M (similar to the construction in [8,Lemmas 3.3,3.4]) and then N( f , u) can be bounded from above by N( f , ϕ). If computed upper bound is small than √ 2 − 1, this implies uniqueness of (u, v).

Lemma 3.3 ([4, Theorem 7]). Let Ω be a bounded polygonal domain in
The remaining parts of this section return to (1.2) with f ∈ L 2 (Ω) and a polygonal domain Ω.

Morley finite element approximation
The first subsection discusses some preliminaries on the Morley FEM and interpolation and enrichment operators. The second subsection derives the existence, uniqueness under a computable smallness assumption and an a priori bound of the discrete solution.

Preliminaries
Let T be an admissible and regular triangulation of the polygonal bounded Lipschitz domain Ω into triangles in R 2 , let h T be the diameter of a triangle T ∈ T and h max := max T ∈T h T . For any ǫ > 0, let T(ǫ) denote the set of all triangulations T with h max < ǫ. For a non-negative integer m, let P m (T ) denote the space of piecewise polynomials of degree at most m. Let Π 0 denote the L 2 projection onto the space P 0 (T ) of piecewise constants and let E and V be the set of edges and vertices of T , respectively. The set of all internal edges (resp. boundary edges) of E is denoted by E(Ω) (resp. E(∂Ω)). Denote the set of internal vertices (resp. vertices on the boundary) of T by V(Ω) (resp. V(∂Ω)). The nonconforming Morley finite element space M(T ) is defined by    [15,19]). There exists a linear operator E M : M(T ) → H 2 0 (Ω) such that any ϕ M ∈ M(T ) satisfies (a)-(d) with a universal constant Λ that depends on the shape-regularity of T but not on the mesh-size h T ∈ P 2 (T ).

Existence, uniqueness, and a priori bound of the discrete solution
This section establishes the well-posedness of the discrete problem (4.1). The discrete analogue This implies The combination of above inequalities, the bound of |||ϕ||| from (2.6) conclude (a) with C d ( . The a priori bound in the equation (4.5) and u M being the minimizer of j pw (

A priori error analysis
This section establishes an a priori error estimates of Morley FEM for the von Kármán obstacle problem with small data.

Main result
Recall M( f , χ) and M d ( f , χ) from Theorems 2.3 and 4.6, the Sobolev and Friedrichs (resp. its discrete versions) constants C S and C F (resp. C dS and C dF ) from (2.1) and Theorem 4.4, and β, α, C(β), and C(α) from Theorem 4.4. The following theorem establishes for small data an a priori energy norm error estimates that is quasi-optimal plus linear convergence.

A priori error analysis of a shifted biharmonic obstacle problem
Let (u, v) be a solution to (1.2) with the regularity u, v ∈ C 2 (Ω) ∩ H 2+α (Ω) ∩ H 2 0 (Ω) from Theorem 3.5. The Sobolev embedding H 2+α (Ω) ֒→ W 2,4 (Ω) leads to f := f + [u, v] ∈ L 2 (Ω). The transformed problem seeks u L ∈ K such that Recall the auxiliary problem from [8] which is a continuous problem with discrete obstacle constraints. Let K A := {ξ ∈ H 2 0 (Ω) : ∀p ∈ V, ξ(p) ≥ χ(p)} for the set V of all vertices in the triangulation T . The biharmonic obstacle problem problem with discrete constraints seeks Equivalently, u A is a minimizer for energy functional J T : The solution u L = u to the biharmonic problem (5.1) and the solution u A to the corresponding auxiliary problem (5.3) satisfy the following result.

Proof of the main result
This and (5.5) imply as h max → 0. Hence, there exists a positive ǫ 2 such that µ d : The later steps of the proof focuses on the error estimates for triangulations T ∈ T(ǫ). Set e := u − u M , δ := v − v M , and the best approximation error RHS := |||u − I M u||| pw + |||v − I M v||| pw + h max .
Step 6 of the proof combines all the estimates T 3 , . . . , T 15 of (5.7) with (5.6). Set positive constants τ 1 := C(1 − µ e (γ + 2)) and τ 2 := 2C(1 − µ e γ −1 ) and deduce with some universal constant C ≈ 1 from the estimates of T 3 , . . . , T 15 . The Young inequality for the last term on the right-hand side of the above estimate imply the assertion This concludes the proof.

Implementation Procedure and Numerical Results
The first subsection is devoted to the implementation procedure to solve the discrete problem (4.1). Subsections 6.2 and 6.3 deal with the results of the numerical experiments and is followed by a subsection on conclusions.
(b) Newtons iteration with initial guess S 0 = (α 0 2 , λ 0 1 , 0, 0). ⋆ For S n := (α n 2 , λ n 1 , β n 1 , β n 2 ), do S n+1 = S n − ∆S n for the solution ∆S n of the linear system of equations J G (S n )∆S n = G(S n ) with J G is the Jacobian matrix of G until ∆S n l 2 (R 2 N ) is less than a given tolerance. • Update m = m + 1. This primal-dual active strategy iteration procedure terminates when Ac m = Ac m−1 and I m = I m−1 .
The flowchart below (see, Figure 1) demonstrates the combined primal-dual active set and Newton algorithms for T 0 , T 1 , . . . . We observe in the examples of this paper (for small f and χ) that at each iteration of primal dual active set algorithm, the Newtons' method converges in four iterations. In this case, we notice that the error between final level and the previous level of the nodal and edge-oriented values in Euclidean norm of R 2N is less than 10 −9 . Also, the primal dual active set algorithm terminates within three steps.

Example 1 (Coincidence set with non-zero measure).
Let the obstacle be given by χ(x) = 1 − 5|x| 2 + |x| 4 , x ∈ Ω = 0.5(−1, 1) 2 . This example is taken from [7]. The discrete coincidence C 6 and C 7 are displayed in Figure 2. Since ∆ 2 χ = 64 > 0 in this example, it is known from [10, Section 8] that the non-coincidence set Ω \ C is connected. This behaviour of the non-coincidence set can be seen in Figure 2 for levels 6 and 7. Table 1 shows errors and orders of convergence for the displacement u and the Airy-stress function v. Observe that linear order of convergences are obtained for u and v in the energy norm, and quadratic order of convergence in L ∞ norm. These numerical order of convergence in energy norm clearly matches the expected order of convergence given in Theorem 5.1. Though the theoretical  rate of convergence in L ∞ norm is not analysed, the numerical rates are obtained similar to that in [7] for the biharmonic obstacle problem.

Example 2. (Coincidence set with zero measure)
In this example taken from [7], χ(x) = 1 − 5|x| 2 − |x| 4 , x ∈ Ω = 0.5(−1, 1) 2 with ∆ 2 χ = −64 < 0 in Ω, and hence, the interior of the coincidence set must be empty, since ∆ 2 u (in the sense of distributions) is a nonnegative measure ([10, Section 8]). This can be observed in the pictures of the discrete coincidence sets displayed in   The errors and orders of convergence for the displacement and the Airy-stress function are presented in Table 2. The orders of convergence results are similar to those obtained in Example 1. Note that Examples 1 and 2 are similar except in the sign of the term |x| 4 that appears in the obstacle function.

The von Kármán obstacle problem on the L-shaped domain
Consider L-shaped domain Ω = (−0.5, 0.5) 2 \ [0, 0.5] 2 , f = 0 and χ(x) = 1 − (x + 0.25) 2 0.2 2 − y 2 0.35 2 as in [7]. Choosing the initial mesh size as h = 0.7071, the successive red-refinement algorithm computes T 1 , . . . , T 5 .  Since Ω is non-convex (reduced elliptic regularity α = 0.5445, [7, Example 4]), we expect only suboptimal order of convergences in energy norm and L ∞ norm, that is, O(h α ) convergence rate in the energy norm (see, Theorem 5.1). However, linear order of convergence is preserved in the energy norm which indicates that the numerical performance is carried out in the non-asymptotic region. The discrete coincidence sets for last two levels are depicted in Figure 4. The non-coincidence set is connected, which agrees with the result in [10] since ∆ 2 χ = 0 in Ω in this example. The convergence rates in Table 3 are not in direct contradiction to Theorem 5.1 but the reduced elliptic regularity suggests a lower rate α = 0.5445 for L-shaped domain. A similar observation is in [7, Table 5.5] with orders of convergence ≈ 0.8 (resp. 1) for energy (resp. L ∞ ) norm. In [7], the numbers are computed with the alternative definitions for error e ℓ (u) := |||u ℓ−1 − u ℓ ||| pw (resp. e ℓ (u) := max p ∈V ℓ−1 |u ℓ−1 (p) − u ℓ (p)|) and order of convergence EOC(ℓ) := log e ℓ−1 (u)/e ℓ (u) /log(2) (resp. log e ℓ−1 (u)/ e ℓ (u) /log (2)). With these definitions, undisplayed numerical experiments confirm the numbers displayed in [7, Table 5.5] precise up to the last digit. This numerical experiment suggests that our implementation is at least consistent with the one in [7]. One possible explanation is that the corner singularity affects the asymptotic convergence rate for very small mesh-sizes only. This is known, for instance, for the L-shaped domain and the Poisson model problem with constant right hand side in the Courant (P 1 conforming) finite element method. The expected rate 2/3 is visible only beyond 2 × 10 6 triangles with far better empirical convergence rates before that.

Conclusions
The numerical results for the Morley FEM in the von Kármán obstacle problem are presented for square domain and L-shaped domain in Sections 6.2 and 6.3. The outputs obtained for the square domain confirm the theoretical rates of convergence given in Theorem 5.1 for α = 1. Example 3 in Section 6.3 illustrates the requirement of smallness assumption on the obstacle for optimal convergence rate. For the L-shaped domain, we expect reduced convergence rates in energy and L ∞ norms from the elliptic regularity. However, linear order of convergence is preserved in the energy norm which indicates that the numerical performance is carried out in the non-asymptotic region.