New High Accuracy Analysis of a Double Set Parameter Nonconforming Element for the Clamped Kirchhoff Plate Unilaterally Constrained by an Elastic Obstacle

: In this paper, a non- C 0 double set parameter ﬁnite element method is presented for the clamped Kirchhoff plate with an elastic unilateral obstacle. A new high accuracy error estimate with order O ( h 2 ) in the broken energy norm is derived by use of a series of novel approaches, including some special features of the element and an incomplete biquadratic interpolation operator. At the same time, some experimental results are provided to verify the theoretical analysis.

It was shown in [1] that the displacement function u can also be derived by the minimization problem: with the total energy J(v) = 1 2 a(v, v) ). Furthermore, the problem (2) is equivalent to the following fourth-order variational inequality of the second kind: In particular, in the limit κ → ∞, the obstacle becomes rigid, and the problem reduces to the constrained minimization:    Find u ∈ K * , such that with K * = {v ∈ V : v ≥ ψ in Ω}, which is the classical displacement obstacle problem of the clamped Kirchhoff plate (cf. [2]) and is equivalent to the following fourth-order variational inequality of the first kind: It is well known that nonconforming finite element methods (FEMs) are an attractive option for solving high order differential equations since the smoothness requirement on finite element functions is weakened. As we know, nonconforming FEMs for the problem (4), the plate obstacle problem with a rigid obstacle, have achieved fruitful results [3][4][5][6][7][8][9], but very few for the elastic obstacle problem (2), except [10]. Compared with the conforming FEMs studied in [1,11], C 0 continuous and discontinuous nonconforming (i.e., C 0 and non-C 0 nonconforming) FEMs were discussed in [10], and convergence analyses were established. On the other hand, the high accuracy analysis of nonconforming FEMs has been an active research area in practical computations, and many high accuracy results have been derived for variational inequalities and the boundary value problem [12][13][14][15][16][17][18][19][20]. However, for the plate obstacle problem (4) with a rigid obstacle, the exact solution u only belongs to H 3 (Ω) ∩ C 2 (Ω) instead of H 4 loc (Ω) [21,22]; thus, the lack of H 4 regularity makes it impossible to develop high accuracy analysis. Fortunately, it was shown in [10] that the solution is more regular if the obstacle is elastic. In fact, the solution of the problem (2) belongs to H 4 (Ω) when the largest interior angles of the domain Ω are smaller than 126.283696... • , and thus, the high accuracy analysis of FEMs for the problem (2) is possible and worth exploring.
In this paper, we attempt to present a high accuracy analysis of nonconforming double set parameter FEMs for the obstacle problem (2). The double set parameter method presented in [23,24] is a useful method to construct unconventional elements, and it involves two sets of parameters. The first set is chosen to meet the convergence requirements according to the generalized patch test [25] or F-E-M-test [26], while the second set is selected to be simple so that the total number of unknowns in the resulting discrete system is small. Up to now, several nonconforming double set parameter plate elements have been successfully applied to deal with the two-sided displacement obstacle problem of the clamped plate [8], plate bending problems [16,27], the linear elasticity problem [28], the fourth-order elliptic singular perturbation problem [29][30][31] and so on. In this paper, a non-C 0 nonconforming double set parameter plate element is employed for the elastic obstacle problem (2). We develop a series of novel approaches including some special features of the element and an incomplete biquadratic interpolation operator so as to get the high accuracy result with order O(h 2 ) in the broken energy norm. Furthermore, we carry out a numerical experiment to show the performance of the proposed method.

A Twelve Parameter Double Set Parameter Element and Its Typical Properties
In the beginning, we introduce a twelve parameter double set parameter element briefly. The readers are referred to [27] for details.
Assume that T h is a regular rectangular decomposition of Ω. For a given K ∈ T h , let its center be (x K , y K ), four vertices be a i (x i , y i ), four edges be F i = a i a i+1 (i = 1, 2, 3, 4(mod4)), and the edge length be 2h x and 2h y , respectively. Moreover, we denote h = max K∈T h {h x , h y }. Besides, assume thatK = [−1, 1] × [−1, 1] is the reference element on the ξ − η plane with a center point (0, 0), four verticesâ 1 = (−1, −1),â 2 = (1, −1),â 3 = (1, 1), andâ 4 = (−1, 1), and four edgesF 1 =â 1â2 ,F 2 =â 2â3 , F 3 =â 3â4 , andF 4 =â 4â1 . Then, there exists an affine mapping J :K −→ K: The shape function space on the element K is taken as: where and P m (K) denotes the space consisting of piecewise polynomials of degree m on element K. The degrees of freedom are selected as: where: For v ∈ P(K), let: where b = (β 1 , β 2 , · · ·, β 12 ) ∈ R 12 , P = (p 1 , p 2 , · · ·, p 12 ) . Substituting (9) into (8), we get: with: It is easy to check that |M| = 2 24 Next, we take nodal parameters as: where v i , v ix , v iy are the function values of v and its first derivatives at vertices a i , respectively. Then, approximating D(v) by a linear combination of the nodal parameters Q(v) as follows: 6,7,8), the numerical integrating values of cubic Hermite interpolation polynomials on the corresponding sides result in: for d i (v) (i = 9, 10, 11, 12), the trapezoidal rule of numerical integration giving: The above discretizations can be rewritten in matrix form as: where: and: (18), it is obvious that the term E(v) does not affect the convergence properties, and we neglect it and introduce a new set of parameters: such that: where , andv x (a) andv y (a) have perturbations on v x (a) and v y (a), respectively. The readers are referred to [23,24] for details. Then, using (10) and (19), we obtain: which in conjunction with (9) leads to the real shape function v. Thus, the unknowns of the discrete system areQ(v), and the real shape function spaceP(K) is a subspace of P(K): It is easy to check that rank(G) = 10, i.e., the dimension ofP(K) is 10, and v ∈P(K) can be expressed Obviously, v ∈P(K) is a cubic polynomial. In other words,P(K) is a completely cubic polynomial space. We denote the corresponding FE space by V h : (R2) F v h ds is continuous across the element edge F and is zero on F ⊂ Γ. (R3) ∀ p(s) ∈ P 1 (F), i.e., p(s) is a linear polynomial function defined on the element edge F, and F p(s) ∂v h ∂n ds is continuous across the element edge F and is zero on Proof. From (8) and (19), we know that D(v h ) is continuous between elements and: Furthermore, it follows from (22) that: which together with (24) and the definition of d i (i = 1, 2, ..., 12) yields (R1), (R2), and (R3). At the same time, suppose v h h = 0, then for any K ∈ T h , it holds that v h | K ∈ P 1 (K). From (R1) and (R2), we can immediately get v h ≡ 0; thus, (R4) holds true.

High Accuracy Analysis
We consider the discrete approximation form of the variational inequality (3) as: where dxdy. It was shown in [1] that the problem (26) is equivalent to the discrete approximation of the plate problem: which is also equivalent to the minimization problem: where the energy functional Since a h (·, ·) is coercive on V h , we can get the lower bound for J h (v h ) and J h (v h ) → ∞ (as v h h → ∞). Then, using the elementary inequality (t − − s − )(t − s) ≥ 0 and an argument similar to that used in [10], we have: In what follows, we will give error estimates for (26).

Theorem 2.
Assume that u and u h are the solutions of (3) and (26), respectively, and the largest internal angle of Ω is ω < 126.283696... • . Then, we have: Here and later, positive constant C is independent of h and may be different at each appearance.
Proof. Employing the triangle inequality yields: where Π h is the associated interpolation operator on V h . Let w h = Π h u − u h . We have by (27) that: Now, we focus on the estimate of the second term on the right-hand side of (31). Obviously, where M(u) = u − (1 − ν) ∂ 2 u ∂s 2 and = (1 − ν). In order to estimate (32), we notice that the solution u of the problem (3) satisfies (see [10]): However, from the construction of the element and (R1)-(R3) in Lemma 1, we know that w h H 1 (Ω).
In this situation, we need to introduce a piecewise interpolation polynomial of w h as follows: Then, from (R1) and (R2) in Lemma 1, we have w I h ∈ H 1 0 (Ω), which in conjunction with (32) and (33) yields: Now, we start to estimate (Er) j one by one for j = 1, 2, ..., 6. Firstly, for φ ∈ C 0 (Ω), let P F 1 φ be the piecewise linear Lagrange interpolation on edge F and P F ∂n ds with | F |= F 1ds, then the corresponding remainders are: Applying (R3) in Lemma 1 leads to: which follows from: Here, the Cauchy-Schwarz inequality and interpolation theorem [32,33] are used. Secondly, for K ∈ T h , we have: Employing integration by parts and (34), the first term of (37) can be estimated as: Similarly, it holds that: which together with interpolation theorem yields: Thirdly, it follows from the definition of w I h in (34) that: Then, let P K 0 ∇( u) = 1 |K| K ∇( u)dxdy and | K |= K 1dxdy. With a similar argument as (Er) 1 and (Er) 2 , we obtain: As for (Er) 4 and (Er) 5 , from the Cauchy-Schwarz inequality and interpolation theorem, it holds that: and: Finally, employing the elementary inequality (t − − s − )(t − s) ≥ 0 and Theorem 1, we have: Combining (31), (35), and the above bounds of (Er) 1 -(Er) 6 results in: Therefore, the desired result (29) follows from (30), (46) and the interpolation theorem immediately.
Remark 1. w I h plays a key role in the estimations of (Er) 2 and (Er) 3 in Theorem 2. If using the error estimate method in [10], we can only obtain the convergence result with the order of O(h). In addition, it should be pointed out that the analysis presented herein is also valid for the elastic obstacle problem with j(v) = Ω κ[(v − ψ) + ] 2 dxdy, t + = max{t, 0}, and f ∈ L 2 (Ω), f ≥ 0 a.e. in Ω.

Remark 2.
Compared with the work in [10], we obtain the same high accuracy result O(h 2 ) with a non-C 0 nonconforming element instead of the C 0 nonconforming ACM's rectangular element. Therefore, our work is an extension of [10], and the requirement for the nonconforming plate element's smoothness from C 0 continuity to mean C 1 continuity is further reduced ((R1)-(R3) implies that the element discussed in this paper is of the mean C 1 type). It can be checked that the high accuracy result in Theorem 2 is no longer true for many classical non-C 0 nonconforming plate elements, such as the Morley triangle element, the Veubeke triangle element, and so on.

Numerical Experiment
In this section, we consider the elastic obstacle problem (2) Table 1. Moreover, Figures 1-3 illustrate the numerical solution u h with N = 64, and Figure 4 presents the errors in the logarithm scales. It is obvious that the errors in the energy norm are convergent at order O(h 2 ), which coincides with the theoretical analysis in Theorem 2.

Conflicts of Interest:
The authors declare no conflicts of interest.