Taylor-Hood like finite elements for nearly incompressible strain gradient elasticity problems

Yulei Liao and Pingbing Ming LSEC, Institute of Computational Mathematics and Scientific/Engineering Computing, AMSS, Chinese Academy of Sciences, No. 55, East Road Zhong-Guan-Cun, Beijing 100190, China School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China [email protected], [email protected] and Yun Xu Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China [email protected]

Abstract.

We propose a family of mixed finite elements that are robust for the nearly incompressible strain gradient model, which is a fourth-order singular perturbed elliptic system. The element is similar to [C. Taylor and P. Hood, Comput. & Fluids, 1(1973), 73-100] in the Stokes flow. Using a uniform discrete B-B inequality for the mixed finite element pairs, we show the optimal rate of convergence that is robust in the incompressible limit. By a new regularity result that is uniform in both the materials parameter and the incompressibility, we prove the method converges with $1/2$ order to the solution with strong boundary layer effects. Moreover, we estimate the convergence rate of the numerical solution to the unperturbed second-order elliptic system. Numerical results for both smooth solutions and the solutions with sharp layers confirm the theoretical prediction.

Key words and phrases:

Mixed finite elements; Nearly incompressible strain gradient elasticity; Uniform error estimate

2020 Mathematics Subject Classification:

Primary 65N15, 65N30; Secondary 74K20

The work of Liao and Ming were supported by National Natural Science Foundation of China through Grant No. 11971467. The work of Xu was supported by National Natural Science Foundation of China through Grant No. 11772067.

1. Introduction

The strain gradient models have drawn great attention recently because they capture the size effect of nano-materials for plasticity [FleckHutchinson:1997] as well as for mechanical meta-materials [Khakalo:2020] by incorporating the higher-order strain gradient and the intrinsic material length scale into the constitutive relations. Studies from the perspective of modeling, mechanics and mathematics may date back to 1960s [Koiter:19641, Mindlin:1964, Hlavacek:19692], while large-scale simulations are relatively recent [zervos:20092, Zybell:2012, Rud:2014, Phunpeng:2015]. Different methods such as H² conforming finite element methods [zervos:20091, Fisher:2010], H¹ conforming mixed finite element methods [Aravas:2002, Phunpeng:2015], nonconforming finite element methods [LiMingShi:2017, LiaoM:2019, LiMingWang:2021], discontinuous Galerkin methods [Engel:2002], isogeometric analysis [Klassen:2011, Niiranen:2016], and meshless methods [Askes:2002] have been used to simulate the strain gradient elastic models with different complexity, just to mention a few. One of the numerical difficulties is that the number of the materials parameters is too large [Mindlin:1964], another is that the materials parameters may cause boundary layer or numerical instability when they tend to certain critical values [Engel:2002].

The strain gradient elasticity model proposed by Altan and Aifantis [Altan:1992] seems the simplest one among them because it contains only one material parameter besides the Lamé constants, while it still models the size effect adequately [Aifantis:1999]. This model is described by the following boundary value problem:

(1.1)

\left\{\begin{aligned} (\iota^{2}\Delta-\boldsymbol{I})\left(\mu\Delta\boldsymbol{u}+(\lambda+\mu)\nabla\operatorname{div}\boldsymbol{u}\right)&=\boldsymbol{f}\quad&&\text{in\;}\Omega,\\ \boldsymbol{u}=\partial_{\boldsymbol{n}}\boldsymbol{u}&=\boldsymbol{0}\quad&&\text{on\;}\partial\Omega,\end{aligned}\right.

where $\Omega\subset\mathbb{R}^{2}$ is a smooth domain, $\boldsymbol{u}:\Omega\to\mathbb{R}^{2}$ is the displacement, $\partial_{\boldsymbol{n}}\boldsymbol{u}$ is the normal derivative of $\boldsymbol{u}$ , $\lambda$ and $\mu$ are the Lamé constants, and $\iota$ is the microscopic parameter such that $0<\iota\leq 1$ , which stands for the intrinsic length scale. Besides modeling the strain gradient elasticity, the system (1.1) may also be regarded as a vector analog of the fourth-order singular perturbed problem, which usually models a clamped plate problem [John:1973, Schuss:1976, Semper:1992, Semper:1994, Tai:2001, Brenner:2011], or arises from a fourth-order perturbation of the fully nonlinear Monge-Ampère equation [Brenner:2011b, Feng:2014]. System of the form (1.1) may also come from the linearized model in MEMS [Laurencot:2017].

In the present work, we are interested in (1.1) for the nearly incompressible materials. Such materials are commonly used in industry and a typical example is natural rubber. To the best of our knowledge, the studies on the approximation of incompressible and nearly incompressible strain gradient elasticity have not been sufficiently addressed in the literature, although vast efforts have been devoted to finite element approximation of the incompressible and nearly incompressible elasticity problems; See e.g., [Hermann:1965, Vog:1983, Simo:1985, Babuska:1992, Brenner:1992, BraessMing:2005, Auricchio:2013]. In [FleckHutchinson:1997]*§III. C, the authors studied the incompressible limit of the strain gradient model. Mixed finite elements for the incompressible Fleck-Hutchinson strain gradient model have been designed and tested in [ShuKingFleck:1999]. A finite element method has been tested for the nearly incompressible strain gradient model in [Wei:2006]. A mixed finite element, which approximated the displacement with Bogner-Fox-Schmidt element [BFS:1965] and approximated the pressure with the $9-$ node quadrilateral element, was constructed for the coupled stress model in [Fisher:2011], and bore certain similarities with problem (1.1). Recently, Hu and Tian [Tian:2021] have proposed several robust elements for the two-dimensional strain gradient model in the framework of reduced integration. Unfortunately, none of the above work proved the robustness of the proposed elements rigorously in the incompressible limit.

Following the classical approach dealing with the nearly incompressible elasticity problem [Hermann:1965], we introduce an auxiliary variable “pressure” $p$ and recast (1.1) into a displacement-pressure mixed variational problem, i.e., the so-called $(\boldsymbol{u},p)-$ formulation. We approximate the displacement by augmenting the finite element space in [GuzmanLeykekhmanNeilan:2012] with certain new bubble functions. The original motivation for the bubble functions is to design the stable finite element pair for the Stokes problem [Arnold:1984, Bernardi:1985]. The augmented bubble functions help out in dealing with the extra constraints such as the divergence stability in Stokes problem and the high order consistency error [Tai:2001, GuzmanLeykekhmanNeilan:2012, Zhang:2012]. Such idea has been exploited by one of the authors to design robust finite elements for the strain gradient elasticity model [LiMingShi:2017]. Besides, we employ the standard continuous Lagrangian finite element of one order lower than that for the displacement to approximate the pressure. Such a finite element pair is similar to the Taylor-Hood element in the Stokes flows [Hood:1973] which is $\mathbb{P}_{r}-\mathbb{P}_{r-1}$ scheme and continuous pressure approximation. For both smooth solutions and solutions with strong boundary layer effects, these mixed finite element pairs are robust in the incompressible limit, here the robustness is understood in the sense that the rate of convergence is uniform in both $\iota$ and $\lambda$ . The bubble function spaces in approximating the displacement are defined by certain orthogonal constraints, and the explicit representations of these spaces are desired for the sake of implementation. We achieve this with the aid of the Jacobi polynomial. In addition to perspicuous results in view of analytics, such representation lends itself to the construction of the analytical shape functions for the approximating space of the displacement. Though we focus on the two-dimensional problem, the element may be readily extended to the three-dimensional problem. cf., Remark 3.3.

By standard mixed finite element theory [BBF:2013], a discrete B-B inequality that is uniform in $\iota$ is needed for the well-posedness of the mixed $(\boldsymbol{u},p)-$ discretization problem. This B-B inequality reduces to the remarkable B-B inequality for the Stokes problem when $\iota$ tends to zero. A natural way to prove the discrete B-B inequality is to construct a uniformly stable Fortin operator [Fortin:1977, Winther:2013, MardalTaiWinther:2002], which does not seem easy due to the complication of the constraints. To this end, we construct a quasi-Fortin operator that takes different forms for small $\iota$ as well as large $\iota$ . This quasi-Fortin operator is bounded in a weighted energy norm in the corresponding regimes of $\iota$ . Besides the discrete B-B inequality, another ingredient in proving the robustness is a new regularity result for the solution of (1.1) that is uniform in both $\lambda$ and $\iota$ , which is crucial to prove the convergence rate for the layered solution. The proof combines the method dealing with the nearly incompressible linear elasticity [Vog:1983] and the regularity estimate for the fourth-order singular perturbed problem [Tai:2001, LiMingWang:2021].

The outline of the paper is as follows. In §2, we introduce Altan and Aifantis’ strain gradient model and its mixed variational formulation. We demonstrate the numerical difficulty caused by large $\lambda$ , and prove the uniform regularity estimate for problem (1.1). In §3, we construct a family of nonconforming finite elements, and derive the explicit formulations for the bubble spaces. In §4, we use the nonconforming elements proposed in §3 together with the continuous Lagrangian finite elements to discretize the mixed variational problem and prove the optimal rate of convergence. In the last section, we report the numerical results, which are consistent with the theoretical prediction.

Throughout this paper, the constant $C$ may differ from line to line, while it is independent of the mesh size $h$ , the materials parameter $\iota$ and the Lamé constant $\lambda$ .

2. The Mixed Variational Formulation and Regularity Estimates

First we fix some notations. The space $L^{2}(\Omega)$ of the square-integrable functions defined on a smooth domain $\Omega$ is equipped with the inner product $(\cdot,\cdot)$ and the norm $\|\,\cdot\,\|_{L^{2}(\Omega)}$ , while $L_{0}^{2}(\Omega)$ is the subspace of $L^{2}(\Omega)$ with mean value zero. Let $H^{m}(\Omega)$ be the standard Sobolev space [AdamsFournier:2003] with the norm $\|\,\cdot\,\|_{H^{m}(\Omega)}$ , while $H_{0}^{m}(\Omega)$ is the closure in $H^{m}(\Omega)$ of $C_{0}^{\infty}(\Omega)$ . We may drop $\Omega$ in $\|\,\cdot\,\|_{H^{m}(\Omega)}$ when no confusion may occur. For any vector-valued function $\boldsymbol{v}$ , its gradient is a matrix-valued function with components $(\nabla\boldsymbol{v})_{ij}=\partial_{j}v_{i}$ , and the symmetric part of $\nabla\boldsymbol{v}$ is defined by $\boldsymbol{\epsilon}(\boldsymbol{v})=(\nabla\boldsymbol{v}+[\nabla\boldsymbol{v}]^{T})/2.$ The divergence operator is defined as $\operatorname{div}\boldsymbol{v}=\partial_{1}v_{1}+\partial_{2}v_{2}$ . The Sobolev spaces $[H^{m}(\Omega)]^{2}$ , $[H_{0}^{m}(\Omega)]^{2}$ and $[L^{2}(\Omega)]^{2}$ of a vector-valued function may be defined similarly as their scalar counterpart. This rule equally applies to the inner products and the norms. The double inner product between tensors $\boldsymbol{A}=(A_{ij})_{i,j=1}^{2},\boldsymbol{B}=(B_{ij})_{i,j=1}^{2}$ equals $\boldsymbol{A}:\boldsymbol{B}=\sum_{i,j=1}^{2}A_{ij}B_{ij}$ .

We recast (1.1) into a variational problem: Find $\boldsymbol{u}\in V{:}=[H^{2}_{0}(\Omega)]^{2}$ such that

(2.1)

a(\boldsymbol{u},\boldsymbol{v})=(\boldsymbol{f},\boldsymbol{v})\quad\text{for all\quad}\boldsymbol{v}\in V,

where $a(\boldsymbol{u},\boldsymbol{v}){:}=(\mathbb{C}\boldsymbol{\epsilon}(\boldsymbol{u}),\boldsymbol{\epsilon}(\boldsymbol{v}))+\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v})),$ and the fourth-order tensor $\mathbb{C}$ and the sixth-order tensor $\mathbb{D}$ are defined as

\mathbb{C}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+2\mu\delta_{ik}\delta_{jl}\quad\text{and}\quad\mathbb{D}_{ijklmn}=\lambda\delta_{il}\delta_{jk}\delta_{mn}+2\mu\delta_{il}\delta_{jm}\delta_{kn},

respectively. Here $\delta_{ij}$ is the Kronecker delta function. The strain gradient $\nabla\boldsymbol{\epsilon}(\boldsymbol{v})$ is a third-order tensor that is defined by $(\nabla\boldsymbol{\epsilon}(\boldsymbol{v}))_{ijk}=\partial_{i}(\boldsymbol{\epsilon}(\boldsymbol{v}))_{jk}$ .

We are interested in the regime when $\lambda\to\infty$ , which means that the material is nearly incompressible. Proceeding along the same line that leads to [LiMingWang:2021]*Theorem 5, the tensor product of the element (NTW) proposed in [Tai:2001] may be used to approximate (1.1), and the error estimate reads as

\|\boldsymbol{u}-\boldsymbol{u}_{h}\|\leq C\lambda(h^{2}+\iota h)\|\,\boldsymbol{u}\,\|_{H^{3}},

where $\|\boldsymbol{v}\|^{2}{:}=a(\boldsymbol{v},\boldsymbol{v})$ , and $C$ is independent of the mesh size $h$ , and $\iota$ and $\lambda$ . Therefore, the error bound degenerates when $\lambda$ is large, and the NTW element does not seem a good candidate for the nearly incompressible material. The following numerical example confirms this observation.

Example 2.1.

Let $\Omega=(0,1)^{2}$ , and $\boldsymbol{u}=(u_{1},u_{2})$ with

u_{1}=-\sin^{3}(\pi x)\sin(2\pi y)\sin(\pi y),\quad u_{2}=\sin(2\pi x)\sin(\pi x)\sin^{3}(\pi y).

It is clear that $\operatorname{div}\boldsymbol{u}=0$ , hence the material is completely incompressible. The details of the numerical experiment such as the mesh generation, are the same as those in § 5. The relative error $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|/\|\boldsymbol{u}\|$ in Table 1 shows that the rate of convergence deteriorates when $\lambda$ is large.

Table 1. Relative errors and rate of convergence for NTW

$\iota\backslash h$	1/8	1/16	1/32	1/64
$\nu=0.3000,\lambda=0.5769,\mu=0.3846$
1e+00	2.681e-01	1.373e-01	6.698e-02	3.334e-02
rate		0.97	1.04	1.01
1e-06	4.550e-02	1.244e-02	3.001e-03	7.467e-04
rate		1.87	2.05	2.01
$\nu=0.499999,\lambda=\text{1.6667e5},\mu=0.3333$
1e+00	9.995e-01	9.979e-01	9.916e-01	9.682e-01
rate		0.00	0.01	0.03
1e-06	9.752e-01	7.233e-01	2.502e-01	6.561e-02
rate		0.43	1.53	1.93

2.1. The mixed variational formulation

We introduce an auxiliary variable $p=\lambda\operatorname{div}\boldsymbol{u},$ and $p\in P{:}=L_{0}^{2}(\Omega)\cap H_{0}^{1}(\Omega)$ . We write Problem (2.1) into a mixed variational problem as

(2.2)

\left\{\begin{aligned} a_{\iota}(\boldsymbol{u},\boldsymbol{v})+b_{\iota}(\boldsymbol{v},p)&=(\boldsymbol{f},\boldsymbol{v})\qquad&&\text{for all\quad}\boldsymbol{v}\in V,\\ b_{\iota}(\boldsymbol{u},q)-\lambda^{-1}c_{\iota}(p,q)&=0\qquad&&\text{for all\quad}q\in P,\end{aligned}\right.

where

$\displaystyle a_{\iota}(\boldsymbol{v},\boldsymbol{w}){:}$	$\displaystyle=2\mu\Bigl{(}(\boldsymbol{\epsilon}(\boldsymbol{v}),\boldsymbol{\epsilon}(\boldsymbol{w}))+\iota^{2}(\nabla\boldsymbol{\epsilon}(\boldsymbol{v}),\nabla\boldsymbol{\epsilon}(\boldsymbol{w}))\Bigr{)},\qquad$	$\displaystyle\boldsymbol{v},\boldsymbol{w}\in V,$
$\displaystyle b_{\iota}(\boldsymbol{v},q){:}$	$\displaystyle=(\operatorname{div}\boldsymbol{v},q)+\iota^{2}(\nabla\operatorname{div}\boldsymbol{v},\nabla q),\qquad$	$\displaystyle\boldsymbol{v}\in V,q\in P,$
$\displaystyle c_{\iota}(s,q){:}$	$\displaystyle=(s,q)+\iota^{2}(\nabla s,\nabla q),\qquad$	$\displaystyle s,q\in P.$

It is convenient to define the weighted norm for all $q\in P$ as

\|q\|_{\iota}{:}=\|\,q\,\|_{L^{2}}+\iota\|\,\nabla q\,\|_{L^{2}}.

$\|q\|_{\iota}$ is a norm over $P$ for any $q\in P$ and any finite $\iota$ . By Poincaré’s inequality, $\|\nabla\boldsymbol{v}\|_{\iota}$ is a norm over $V$ for any $\boldsymbol{v}\in V$ . To study the well-posedness of Problem (2.2), we start with the following B-B inequality that is uniform for any $\iota$ .

Lemma 2.2.

For any $q\in P$ , there exists $\boldsymbol{v}\in V$ such that

(2.3)

\operatorname{div}\boldsymbol{v}=q\qquad\text{and}\qquad\|\nabla\boldsymbol{v}\|_{\iota}\leq C\|q\|_{\iota},

where $C$ only depends on $\Omega$ but is independent of $\iota$ .

Proof.

By [Galdi:2011]*Theorem III 3.3 and [Costabel:2010]*Proposition 4.1, for any $q\in P$ , there exists $\boldsymbol{v}\in V$ such that $\operatorname{div}\boldsymbol{v}=q$ and

(2.4)

\|\,\boldsymbol{v}\,\|_{H^{1}}\leq C\|\,q\,\|_{L^{2}}\qquad\text{and}\qquad\|\,\boldsymbol{v}\,\|_{H^{2}}\leq C\|\,q\,\|_{H^{1}},

where the constant $C$ only depends on $\Omega$ .

Because the mean of $q$ is zero for any $q\in P$ , by Poincaré’s inequality, there exists $C$ such that

\|\,q\,\|_{H^{1}}\leq C\|\,\nabla q\,\|_{L^{2}}.

Combining the above two inequalities, we obtain

\|\nabla\boldsymbol{v}\|_{\iota}=\|\,\nabla\boldsymbol{v}\,\|_{L^{2}}+\iota\|\,\nabla^{2}\boldsymbol{v}\,\|_{L^{2}}\leq C\|q\|_{\iota}.

This gives (2.3). ∎

Lemma 2.3.

There exists a unique $\boldsymbol{u}\in V$ and $p\in P$ satisfying (2.2), and there exists $C$ independent of $\iota$ and $\lambda$ such that

(2.5)

\|\nabla\boldsymbol{u}\|_{\iota}+\|p\|_{\iota}\leq C\|\,\boldsymbol{f}\,\|_{H^{-1}}.

Proof.

By the first Korn’s inequality [Korn:1908, Korn:1909],

(2.6)

\|\,\boldsymbol{\epsilon}(\boldsymbol{v})\,\|_{L^{2}}^{2}\geq\dfrac{1}{2}\|\,\nabla\boldsymbol{v}\,\|_{L^{2}}^{2}\qquad\text{for all\quad}\boldsymbol{v}\in[H_{0}^{1}(\Omega)]^{2},

and the H² Korn’s inequality [LiMingWang:2021]*Theorem 1,

(2.7)

\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{v})\,\|_{L^{2}}^{2}\geq\left(1-\dfrac{1}{\sqrt{2}}\right)\|\,\nabla^{2}\boldsymbol{v}\,\|_{L^{2}}^{2}\qquad\text{for all\quad}\boldsymbol{v}\in[H^{2}(\Omega)]^{2},

we obtain

a_{\iota}(\boldsymbol{v},\boldsymbol{v})\geq 2\mu\left(\dfrac{1}{2}\|\,\nabla\boldsymbol{v}\,\|_{L^{2}}^{2}+\Bigl{(}1-\dfrac{1}{\sqrt{2}}\Bigr{)}\iota^{2}\|\,\nabla^{2}\boldsymbol{v}\,\|_{L^{2}}^{2}\right)\geq\dfrac{\mu}{2}\|\nabla\boldsymbol{v}\|_{\iota}^{2}.

Using (2.3), for any $p\in P$ , there exists $\boldsymbol{v}_{0}\in V$ such that $\operatorname{div}\boldsymbol{v}_{0}=p$ and $\|\nabla\boldsymbol{v}_{0}\|_{\iota}\leq C\|p\|_{\iota}$ . This implies

\sup_{\boldsymbol{v}\in V}\dfrac{b_{\iota}(\boldsymbol{v},p)}{\|\nabla\boldsymbol{v}\|_{\iota}}\geq\dfrac{b_{\iota}(\boldsymbol{v}_{0},p)}{\|\nabla\boldsymbol{v}_{0}\|_{\iota}}=\dfrac{\|p\|_{\iota}^{2}}{\|\nabla\boldsymbol{v}_{0}\|_{\iota}}\geq C\|p\|_{\iota}.

By [Braess:1996]*Theorem 2, we immediately obtain the well-posedness of (2.2) and the estimate (2.5) by noting

\lvert(\boldsymbol{f},\boldsymbol{v})\rvert\leq\|\,\boldsymbol{f}\,\|_{H^{-1}}\|\,\boldsymbol{v}\,\|_{H^{1}}\leq C\|\,\boldsymbol{f}\,\|_{H^{-1}}\|\,\nabla\boldsymbol{v}\,\|_{L^{2}}\leq C\|\,\boldsymbol{f}\,\|_{H^{-1}}\|\nabla\boldsymbol{v}\|_{\iota}.

∎

By the standard regularity theory for the elliptic system, we find $\boldsymbol{u}\in[H^{4}(\Omega)]^{2}$ and $p\in H^{3}(\Omega)$ provided that $\boldsymbol{f}\in[L^{2}(\Omega)]^{2}$ , while we are interested in whether the shift estimate $\|\,\nabla^{2}\boldsymbol{u}\,\|_{\iota}+\|\,\nabla p\,\|_{\iota}\leq C(\iota)\|\,\boldsymbol{f}\,\|_{L^{2}}$ is true with a $\lambda-$ independent constant $C(\iota)$ , this is the objective of the next part.

2.2. Regularity estimates

We aim to study the regularity of the solution of (1.1). Letting $\iota\to 0$ , we find $\boldsymbol{u}_{0}\in[H_{0}^{1}(\Omega)]^{2}$ satisfying

(2.8)

-\mathcal{L}\boldsymbol{u}_{0}=\boldsymbol{f}\quad\text{in\;}\Omega,\qquad\boldsymbol{u}_{0}=\boldsymbol{0}\quad\text{on\;}\partial\Omega,

in the sense of distribution, where $\mathcal{L}\boldsymbol{u}_{0}{:}=\mu\Delta\boldsymbol{u}_{0}+(\lambda+\mu)\nabla\operatorname{div}\boldsymbol{u}_{0}$ . The H ${}^{1}-$ error for $\boldsymbol{u}-\boldsymbol{u}_{0}$ will be given in Theorem 2.9, which is crucial for the regularity estimate of problem (1.1). We reshape (2.8) into a variational problem: Find $\boldsymbol{u}_{0}\in[H_{0}^{1}(\Omega)]^{2}$ such that

(2.9)

(\mathbb{C}\boldsymbol{\epsilon}(\boldsymbol{u}_{0}),\boldsymbol{\epsilon}(\boldsymbol{v}))=(\boldsymbol{f},\boldsymbol{v})\qquad\text{for all\;}\boldsymbol{v}\in[H_{0}^{1}(\Omega)]^{2}.

By [BacutaBramble:2003], we have the following shift estimate for $\boldsymbol{u}_{0}$ : There exists $C$ independent of $\lambda$ such that

(2.10)

\|\,\boldsymbol{u}_{0}\,\|_{H^{2}}+\lambda\|\,\operatorname{div}\boldsymbol{u}_{0}\,\|_{H^{1}}\leq C\|\,\boldsymbol{f}\,\|_{L^{2}}.

Next we study an auxiliary boundary value problem:

(2.11)

\left\{\begin{aligned} \Delta\mathcal{L}\boldsymbol{w}&=\boldsymbol{F},\qquad\text{in\quad}\Omega,\\ \boldsymbol{w}=\partial_{\boldsymbol{n}}\boldsymbol{w}&=\boldsymbol{0},\qquad\text{on\quad}\partial\Omega.\end{aligned}\right.

The a-priori estimate for the solution of the above boundary value problem reads as

Lemma 2.4.

There exists a unique $\boldsymbol{w}\in V$ satisfying (2.11), and there exists $C$ independent of $\lambda$ such that

(2.12)

\|\,\boldsymbol{w}\,\|_{H^{2}}+\lambda\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{1}}\leq C\|\,\boldsymbol{F}\,\|_{H^{-2}}.

Proof.

We recast (2.11) into a variational problem: Find $\boldsymbol{w}\in V$ such that

A(\boldsymbol{w},\boldsymbol{v})=(\boldsymbol{F},\boldsymbol{v})\qquad\text{for all\quad}\boldsymbol{v}\in V,

where $A(\boldsymbol{v},\boldsymbol{z}){:}=2\mu(\nabla\boldsymbol{\epsilon}(\boldsymbol{v}),\nabla\boldsymbol{\epsilon}(\boldsymbol{z}))+\lambda(\nabla\operatorname{div}\boldsymbol{v},\nabla\operatorname{div}\boldsymbol{z})$ for any $\boldsymbol{v},\boldsymbol{z}\in V$ .

For any $\boldsymbol{v}\in V$ , by the H ${}^{2}-$ Korn’s inequality (2.7) and Poincaré’s inequality, there exists $C$ such that

A(\boldsymbol{v},\boldsymbol{v})\geq 2\mu\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{v})\,\|_{L^{2}}^{2}\geq\dfrac{\mu}{2}\|\,\nabla^{2}\boldsymbol{v}\,\|_{L^{2}}^{2}\geq C\|\,\boldsymbol{v}\,\|_{H^{2}}^{2}.

The existence and uniqueness of $\boldsymbol{w}\in V$ follow from the Lax-Milgram theorem, and

(2.13)

\|\,\nabla^{2}\boldsymbol{w}\,\|_{L^{2}}\leq\|\,\boldsymbol{w}\,\|_{H^{2}}\leq C\|\,\boldsymbol{F}\,\|_{H^{-2}}.

Noting that $\operatorname{div}\boldsymbol{w}\in P$ , using (2.4), we obtain that, there exists $\boldsymbol{v}_{0}\in V$ such that $\operatorname{div}\boldsymbol{v}_{0}=\operatorname{div}\boldsymbol{w},$ and

\|\,\nabla^{2}\boldsymbol{v}_{0}\,\|_{L^{2}}\leq C\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{1}}\leq C\|\,\nabla\operatorname{div}\boldsymbol{w}\,\|_{L^{2}}.

A combination of the above two inequalities gives

	$\displaystyle\lambda\\|\,\nabla\operatorname{div}\boldsymbol{w}\,\\|_{L^{2}}^{2}$	$\displaystyle=\lambda(\nabla\operatorname{div}\boldsymbol{w},\nabla\operatorname{div}\boldsymbol{v}_{0})=A(\boldsymbol{w},\boldsymbol{v}_{0})-2\mu(\nabla\boldsymbol{\epsilon}(\boldsymbol{w}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))$
		$\displaystyle=(\boldsymbol{F},\boldsymbol{v}_{0})-2\mu(\nabla\boldsymbol{\epsilon}(\boldsymbol{w}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))$
		$\displaystyle\leq\\|\,\boldsymbol{F}\,\\|_{H^{-2}}\\|\,\boldsymbol{v}_{0}\,\\|_{H^{2}}+2\mu\\|\,\nabla^{2}\boldsymbol{w}\,\\|_{L^{2}}\\|\,\nabla^{2}\boldsymbol{v}_{0}\,\\|_{L^{2}}$
		$\displaystyle\leq C\left(\\|\,\boldsymbol{F}\,\\|_{H^{-2}}+2\mu\\|\,\nabla^{2}\boldsymbol{w}\,\\|_{L^{2}}\right)\\|\,\nabla^{2}\boldsymbol{v}_{0}\,\\|_{L^{2}}$
		$\displaystyle\leq C\\|\,\boldsymbol{F}\,\\|_{H^{-2}}\\|\,\nabla\operatorname{div}\boldsymbol{w}\,\\|_{L^{2}}.$

This implies $\lambda\|\,\nabla\operatorname{div}\boldsymbol{w}\,\|_{L^{2}}\leq C\|\,\boldsymbol{F}\,\|_{H^{-2}},$ which together with (2.13) and Poincaré’s inequality gives (2.12). ∎

Now we turn to prove the regularity estimate of problem (2.11). We consider an auxiliary elliptic system: For any $\widetilde{\boldsymbol{F}}\in[L^{2}(\Omega)]^{2}$ and $\widetilde{G}\in H^{1}(\Omega)$ , find $\boldsymbol{z}\in V$ and $q\in P$ such that the following boundary value problem is valid in the sense of distribution:

(2.14)

\left\{\begin{aligned} \mu\Delta^{2}\boldsymbol{z}+\nabla\Delta q&=\widetilde{\boldsymbol{F}}\qquad&&\text{in\quad}\Omega,\\ \Delta\operatorname{div}\boldsymbol{z}&=\widetilde{G}\qquad&&\text{in\quad}\Omega,\\ \boldsymbol{z}=\partial_{\boldsymbol{n}}\boldsymbol{z}&=\boldsymbol{0}\qquad&&\text{on\quad}\partial\Omega,\\ q&=0\qquad&&\text{on\quad}\partial\Omega.\end{aligned}\right.

Lemma 2.5.

Let $\boldsymbol{z}\in[H^{4}(\Omega)]^{2}$ and $q\in H^{3}(\Omega)$ be the solution of (2.14). Assume that $m$ is a nonnegative integer, then there exists $C$ depending only on $\Omega$ and $\mu$ such that

(2.15)

\|\,\boldsymbol{z}\,\|_{H^{m+4}}+\|\,q\,\|_{H^{m+3}}\leq C\left(\|\,\widetilde{\boldsymbol{F}}\,\|_{H^{m}}+\|\,\widetilde{G}\,\|_{H^{m+1}}+\|\,\boldsymbol{z}\,\|_{L^{2}}+\|\,q\,\|_{L^{2}}\right).

Proof.

We write (2.14)₁ and (2.14)₂ as

\begin{pmatrix}\mu\Delta^{2}&0&\partial_{x}\Delta\\ 0&\mu\Delta^{2}&\partial_{y}\Delta\\ \Delta\partial_{x}&\Delta\partial_{y}&0\end{pmatrix}\begin{pmatrix}z_{1}\\ z_{2}\\ q\end{pmatrix}=\begin{pmatrix}\widetilde{F}_{1}\\ \widetilde{F}_{2}\\ \widetilde{G}\end{pmatrix}.

The symbol of the above system is

\mathcal{L}(\xi)=\begin{pmatrix}\mu\lvert\xi\rvert^{4}&0&\xi_{1}\lvert\xi\rvert^{2}\\ 0&\mu\lvert\xi\rvert^{4}&\xi_{2}\lvert\xi\rvert^{2}\\ \xi_{1}\lvert\xi\rvert^{2}&\xi_{2}\lvert\xi\rvert^{2}&0\end{pmatrix}.

A direct calculation gives

\lvert\det\mathcal{L}(\xi)\rvert=\mu\lvert\xi\rvert^{10}>0\qquad\text{if\quad}\xi\not=0.

This means that the boundary value problem (2.14) is elliptic in the sense of Agmon-Douglis-Nirenberg [Agmon:1964]. Moreover, the boundary condition is pure Dirichlet, and it is straightforward to verify that the boundary condition satisfies the complementing condition [Agmon:1964]. The regularity estimate (2.15) follows from [Agmon:1964]. ∎

A direct consequence of the above lemma is the following regularity estimate for problem (2.11).

Lemma 2.6.

Let $\boldsymbol{w}\in V$ be the solution of (2.11), there exists $C$ independent of $\lambda$ such that

(2.16)

\|\,\boldsymbol{w}\,\|_{H^{3}}+\lambda\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{2}}\leq C\|\,\boldsymbol{F}\,\|_{H^{-1}}.

Proof.

Using the standard elliptic regularity estimate, there exists a unique solution $\boldsymbol{w}\in[H^{4}(\Omega)]^{2}$ when $\boldsymbol{F}\in[L^{2}(\Omega)]^{2}$ . We introduce $\boldsymbol{z}=\boldsymbol{w}$ and $q=(\lambda+\mu)\operatorname{div}\boldsymbol{w}$ , hence $\boldsymbol{z}\in[H^{4}(\Omega)]^{2}$ and $q\in H^{3}(\Omega)$ satisfy (2.14) with $\widetilde{\boldsymbol{F}}=\boldsymbol{F}$ and $\widetilde{G}=\Delta\operatorname{div}\boldsymbol{w}$ .

By (2.15) with $m=0$ , we obtain

\|\,\boldsymbol{w}\,\|_{H^{4}}+(\lambda+\mu)\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}\leq C\left(\|\,\boldsymbol{F}\,\|_{L^{2}}+\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}+\|\,\boldsymbol{w}\,\|_{L^{2}}+(\lambda+\mu)\|\,\operatorname{div}\boldsymbol{w}\,\|_{L^{2}}\right).

Using the a-priori estimate (2.12), we obtain

\|\,\boldsymbol{w}\,\|_{H^{4}}+(\lambda+\mu)\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}\leq C_{0}\left(\|\,\boldsymbol{F}\,\|_{L^{2}}+\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}\right).

Now for $\lambda+\mu>2C_{0}$ , it follows from

\|\,\boldsymbol{w}\,\|_{H^{4}}+(\lambda+\mu)\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}\leq C_{0}\|\,\boldsymbol{F}\,\|_{L^{2}}+\dfrac{\lambda+\mu}{2}\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}

that

\|\,\boldsymbol{w}\,\|_{H^{4}}+\lambda\|\,\operatorname{div}\boldsymbol{w}\,\|_{H^{3}}\leq 2C_{0}\|\,\boldsymbol{F}\,\|_{L^{2}}.

Interpolating the above inequality with (2.12), we obtain (2.16).

If $\lambda+\mu\leq 2C_{0}$ , then (2.16) follows from the standard regularity estimates [Agmon:1964] for problem (2.11). ∎

We turn to prove the regularity of problem (1.1) when $\boldsymbol{f}\in[L^{2}(\Omega)]^{2}$ . Let $\boldsymbol{u}$ and $\boldsymbol{u}_{0}$ be the solutions of (1.1) and (2.8), respectively. For any $\boldsymbol{v}\in[H_{0}^{1}(\Omega)\cap H^{2}(\Omega)]^{2}$ , integration by parts gives

(2.17)

(\mathbb{C}\boldsymbol{\epsilon}(\boldsymbol{u}),\boldsymbol{\epsilon}(\boldsymbol{v}))=-(\mathcal{L}\boldsymbol{u},\boldsymbol{v}),

and

(2.18)

\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}))=\iota^{2}(\Delta\mathcal{L}\boldsymbol{u},\boldsymbol{v})+\iota^{2}\int_{\partial\Omega}{(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})}\cdot\partial_{\boldsymbol{n}}\boldsymbol{v}\mathrm{d}\sigma(\boldsymbol{x}),

where $\boldsymbol{\sigma}{:}=2\mu\boldsymbol{\epsilon}(\boldsymbol{u})+\lambda\operatorname{div}\boldsymbol{u}\boldsymbol{I}$ and $(\partial_{\boldsymbol{n}}\boldsymbol{\sigma})_{ij}{:}=\partial_{\boldsymbol{n}}\sigma_{ij}$ . The boundary term in (2.18) is derived by the fact $\partial_{j}v_{i}=n_{j}\partial_{\boldsymbol{n}}v_{i}+t_{j}\partial_{\boldsymbol{t}}v_{i}=n_{j}\partial_{\boldsymbol{n}}v_{i}$ and

(2.19)		$\displaystyle 2\mu\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v})\mathrm{d}\sigma(\boldsymbol{x})$	$\displaystyle+\lambda\int_{\partial\Omega}\partial_{\boldsymbol{n}}\operatorname{div}\boldsymbol{u}\operatorname{div}\boldsymbol{v}\mathrm{d}\sigma(\boldsymbol{x})=\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\sigma}:\nabla\boldsymbol{v}\mathrm{d}\sigma(\boldsymbol{x})$
(2.19)		$\displaystyle=$	$\displaystyle\int_{\partial\Omega}\partial_{\boldsymbol{n}}\sigma_{ij}\partial_{j}v_{i}\mathrm{d}\sigma(\boldsymbol{x})=\int_{\partial\Omega}\partial_{\boldsymbol{n}}\sigma_{ij}n_{j}\partial_{\boldsymbol{n}}v_{i}\mathrm{d}\sigma(\boldsymbol{x})=\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}v\mathrm{d}\sigma(\boldsymbol{x}).$

A combination of (2.17) and (2.18) leads to

(\mathbb{C}\boldsymbol{\epsilon}(\boldsymbol{u}),\boldsymbol{\epsilon}(\boldsymbol{v}))+\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}))=\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{v}\mathrm{d}\sigma(\boldsymbol{x})+(\boldsymbol{f},\boldsymbol{v}),

which together with (2.9) yields

(2.20)

(\mathbb{C}\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}),\boldsymbol{\epsilon}(\boldsymbol{v}))+\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}))=\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{v}\mathrm{d}\sigma(\boldsymbol{x}).

This identity is the starting point of the proof.

We shall frequently use the following multiplicative trace inequality. There exists $C$ that depends only on $\Omega$ such that

(2.21)

\|\,\psi\,\|_{L^{2}(\partial\Omega)}\leq C\|\,\psi\,\|_{L^{2}(\Omega)}^{1/2}\|\,\psi\,\|_{H^{1}(\Omega)}^{1/2},\qquad\psi\in H^{1}(\Omega).

Using Lemma 2.6, we condense the regularity of problem (1.1) to estimating $\boldsymbol{u}-\boldsymbol{u}_{0}$ and $p-p_{0}$ in various norms, with $p=\lambda\operatorname{div}\boldsymbol{u}$ and $p_{0}=\lambda\operatorname{div}\boldsymbol{u}_{0}$ .

Lemma 2.7.

There exists $C$ independent of $\iota$ and $\lambda$ such that

(2.22)

\|\,\boldsymbol{u}\,\|_{H^{3}}+\|\,p\,\|_{H^{2}}\leq C\iota^{-2}(\|\,\epsilon(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}}+\|\,p-p_{0}\,\|_{L^{2}}).

Proof.

We rewrite (1.1) as

\left\{\begin{aligned} \Delta\mathcal{L}\boldsymbol{u}&=\iota^{-2}\mathcal{L}(\boldsymbol{u}-\boldsymbol{u}_{0})\qquad&&\text{in\quad}\Omega,\\ \boldsymbol{u}=\partial_{\boldsymbol{n}}\boldsymbol{u}&=\boldsymbol{0}\qquad&&\text{on\quad}\partial\Omega.\end{aligned}\right.

Applying the regularity estimate (2.16) to the elliptic system (2.11), and using (2.17), we obtain

\|\,\boldsymbol{u}\,\|_{H^{3}}+\|\,p\,\|_{H^{2}}\leq C\iota^{-2}\|\,\mathcal{L}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{H^{-1}}\leq C\iota^{-2}(\|\,\epsilon(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}}+\|\,p-p_{0}\,\|_{L^{2}}).

∎

The next lemma is crucial to prove Theorem 2.9, which transforms the estimate of $\|\,p-p_{0}\,\|_{\iota}$ in terms of $\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}$ besides a term concerning $\boldsymbol{f}$ .

Lemma 2.8.

There exists $C$ independent of $\iota$ and $\lambda$ such that

(2.23)

\|\,p-p_{0}\,\|_{\iota}\leq C(\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}+\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}).

Proof.

By [Danchin:2013]*Theorem 3.1 and Poincaré’s inequality, there exists $\boldsymbol{v}_{0}\in[H^{2}(\Omega)\cap H_{0}^{1}(\Omega)]^{2}$ such that $\operatorname{div}\boldsymbol{v}_{0}=\operatorname{div}(\boldsymbol{u}-\boldsymbol{u}_{0})$ , and

(2.24)

\|\,\boldsymbol{v}_{0}\,\|_{H^{1}}\leq C\|\,\operatorname{div}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}},\qquad\|\,\boldsymbol{v}_{0}\,\|_{H^{2}}\leq C\|\,\nabla\operatorname{div}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}}.

Substituting $\boldsymbol{v}=\boldsymbol{v}_{0}$ into (2.20), multiplying the resulting identity by $\lambda$ , we obtain

	$\displaystyle\\|\,p-p_{0}\,\\|_{\iota}^{2}$	$\displaystyle=-\iota^{2}(\nabla p_{0},\nabla(p-p_{0}))-2\lambda\mu\Bigl{(}(\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}),\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))+\iota^{2}(\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))\Bigr{)}$
(2.25)			$\displaystyle\quad+\lambda\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{v}_{0}\mathrm{d}\sigma(\boldsymbol{x}).$

By the regularity estimate (2.10), the first term may be bounded as

\iota^{2}\lvert(\nabla p_{0},\nabla(p-p_{0}))\rvert\leq\dfrac{\iota^{2}}{8}\|\,\nabla(p-p_{0})\,\|_{L^{2}}^{2}+2\iota^{2}\|\,\nabla p_{0}\,\|_{L^{2}}^{2}\leq\dfrac{1}{8}\|\,p-p_{0}\,\|_{\iota}^{2}+C\iota^{2}\|\,\boldsymbol{f}\,\|_{L^{2}}^{2}.

Using the triangle inequality and (2.10) again, we obtain

(2.26)

\iota\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\|_{L^{2}}\leq\iota\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}}+\iota\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u}_{0})\,\|_{L^{2}}\leq\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}+C\iota\|\,\boldsymbol{f}\,\|_{L^{2}}.

Using (2.24) and the above inequality, we bound the second term as

	$\displaystyle 2\lambda\mu\lvert(\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}),\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))+\iota^{2}(\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))\rvert$	$\displaystyle\leq C(\\|\,\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{L^{2}}+\iota\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}})\\|\,p-p_{0}\,\\|_{\iota}$
		$\displaystyle\leq\dfrac{1}{8}\\|\,p-p_{0}\,\\|_{\iota}^{2}+C\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota^{2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}\right).$

Using (2.19) and the definition of $\boldsymbol{v}_{0}$ , we rewrite the boundary term as

	$\displaystyle\lambda\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{v}_{0}\mathrm{d}\sigma(\boldsymbol{x})$	$\displaystyle=2\lambda\mu\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})+\lambda^{2}\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\operatorname{div}\boldsymbol{u}\operatorname{div}\boldsymbol{v}_{0}\mathrm{d}\sigma(\boldsymbol{x})$
		$\displaystyle=2\lambda\mu\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})+\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}p(p-p_{0})\mathrm{d}\sigma(\boldsymbol{x})$
		$\displaystyle=2\lambda\mu\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})-\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}p\;p_{0}\mathrm{d}\sigma(\boldsymbol{x}).$

Recalling the trace inequality (2.21), we estimate the first boundary term as

	$\displaystyle 2\lambda\mu\iota^{2}\biggl{\|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})\biggr{\|}$	$\displaystyle\leq 2\lambda\mu\iota^{2}\\|\,\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}(\partial\Omega)}\\|\,\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\,\\|_{L^{2}(\partial\Omega)}$
		$\displaystyle\leq C\lambda\iota^{2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{H^{1}}^{1/2}\\|\,\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\,\\|_{L^{2}}^{1/2}\\|\,\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\,\\|_{H^{1}}^{1/2}.$

Using (2.24), there exists $C$ independent of $\lambda$ and $\iota$ such that

\lambda^{2}\|\,\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\,\|_{L^{2}}\|\,\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\,\|_{H^{1}}\leq C\|\,p-p_{0}\,\|_{L^{2}}\|\,\nabla(p-p_{0})\,\|_{L^{2}}\leq C\iota^{-1}\|\,p-p_{0}\,\|_{\iota}^{2}.

Using (2.26) to estimate $\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\|_{L^{2}}$ and using (2.22) to bound $\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\|_{H^{1}}$ , we obtain

	$\displaystyle\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{H^{1}}$	$\displaystyle\leq C\iota^{-3}\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}}\right)\left(\\|\,\epsilon(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{L^{2}}+\\|\,p-p_{0}\,\\|_{L^{2}}\right)$
		$\displaystyle\leq C\iota^{-3}\Bigl{(}\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota^{2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}+\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}}\right)\\|\,p-p_{0}\,\\|_{\iota}\Bigr{)}.$

A combination of the above three inequalities gives

	$\displaystyle 2\lambda\mu\iota^{2}\biggl{\|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})\biggr{\|}\leq$	$\displaystyle C\Bigl{(}(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}})\\|\,p-p_{0}\,\\|_{\iota}$
		$\displaystyle\qquad+(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}})^{1/2}\\|\,p-p_{0}\,\\|_{\iota}^{3/2}\Bigr{)}$
	$\displaystyle\leq$	$\displaystyle\dfrac{1}{8}\\|\,p-p_{0}\,\\|_{\iota}^{2}+C\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota^{2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}\right).$

Using the trace inequality (2.21) and the regularity estimate (2.10) again, we bound

\iota^{2}\biggl{|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}pp_{0}\mathrm{d}\sigma(\boldsymbol{x})\biggr{|}\leq C\iota^{2}\|\,\nabla p\,\|_{L^{2}}^{1/2}\|\,\nabla p\,\|_{H^{1}}^{1/2}\|\,p_{0}\,\|_{H^{1}}\leq C\iota^{2}\|\,\nabla p\,\|_{L^{2}}^{1/2}\|\,\nabla p\,\|_{H^{1}}^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}.

Using the triangle inequality and (2.10) again, we obtain

\iota\|\,\nabla p\,\|_{L^{2}}\leq\|\,p-p_{0}\,\|_{\iota}+\iota\|\,\nabla p_{0}\,\|_{L^{2}}\leq\|\,p-p_{0}\,\|_{\iota}+C\iota\|\,\boldsymbol{f}\,\|_{L^{2}},

which together with (2.22) implies

\|\,\nabla p\,\|_{L^{2}}\|\,\nabla p\,\|_{H^{1}}\leq C\iota^{-3}\Bigl{(}\|\,p-p_{0}\,\|_{\iota}^{2}+(\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}+\iota\|\,\boldsymbol{f}\,\|_{L^{2}})\|\,p-p_{0}\,\|_{\iota}+\iota\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}\|\,\boldsymbol{f}\,\|_{L^{2}}\Bigr{)}.

Combining the above three inequalities, we bound the second boundary term as

	$\displaystyle\iota^{2}\biggl{\|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}pp_{0}\mathrm{d}\sigma(\boldsymbol{x})\biggr{\|}\leq$	$\displaystyle C\Bigl{(}\iota^{1/2}\\|\,p-p_{0}\,\\|_{\iota}\\|\,\boldsymbol{f}\,\\|_{L^{2}}+(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}})^{3/2}\\|\,p-p_{0}\,\\|_{\iota}^{1/2}$
		$\displaystyle\qquad+\iota\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{3/2}\Bigr{)}$
		$\displaystyle\leq\dfrac{1}{8}\\|\,p-p_{0}\,\\|_{\iota}^{2}+C\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}\right).$

Substituting the above inequalities into (2.25), we obtain

\|\,p-p_{0}\,\|_{\iota}^{2}\leq\dfrac{1}{2}\|\,p-p_{0}\,\|_{\iota}^{2}+C\left(\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}^{2}+\iota\|\,\boldsymbol{f}\,\|_{L^{2}}^{2}\right).

This immediately gives (2.23). ∎

We are ready to prove the regularity of problem (1.1).

Theorem 2.9.

There exists $C$ independent of $\iota$ and $\lambda$ such that

(2.27)

\|\,\boldsymbol{u}-\boldsymbol{u}_{0}\,\|_{H^{1}}+\|\,p-p_{0}\,\|_{L^{2}}\leq C\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}},

and

(2.28)

\|\,\boldsymbol{u}\,\|_{H^{2+k}}+\|\,p\,\|_{H^{1+k}}\leq C\iota^{-1/2-k}\|\,\boldsymbol{f}\,\|_{L^{2}},\qquad k=0,1.

The estimate (2.27) improves the known results [LiMingWang:2021]*Lemma 1 in two aspects. It clarifies the fact that the estimate is $\lambda-$ independent and it gives the convergence rate for the pressure. The rate $\iota^{1/2}$ is optimal even for the scalar counterpart; cf., [Tai:2001]*Lemma 5.1.

Proof.

Substituting $\boldsymbol{v}=\boldsymbol{u}-\boldsymbol{u}_{0}$ into (2.20), we get

a(\boldsymbol{u}-\boldsymbol{u}_{0},\boldsymbol{u}-\boldsymbol{u}_{0})=-\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}))-\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{u}_{0}\mathrm{d}\sigma(\boldsymbol{x}).

Using the regularity estimate (2.10), we bound the first term as

\iota^{2}\lvert(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}),\nabla\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}))\rvert\leq\dfrac{\iota^{2}}{4}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}),\nabla\boldsymbol{\epsilon}(\boldsymbol{u}-\boldsymbol{u}_{0}))+\iota^{2}(\mathbb{D}\nabla\boldsymbol{\epsilon}(\boldsymbol{u}_{0}),\nabla\boldsymbol{\epsilon}(\boldsymbol{u}_{0}))\leq\dfrac{1}{4}a(\boldsymbol{u}-\boldsymbol{u}_{0},\boldsymbol{u}-\boldsymbol{u}_{0})+C\iota^{2}\|\,\boldsymbol{f}\,\|_{L^{2}}^{2}.

To bound the second term, we let $\boldsymbol{v}=\boldsymbol{u}_{0}$ in (2.19) and obtain

\iota^{2}\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{u}_{0}\mathrm{d}\sigma(\boldsymbol{x})=2\mu\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{u}_{0})\mathrm{d}\sigma(\boldsymbol{x})+\iota^{2}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\operatorname{div}\boldsymbol{u}\;p_{0}\mathrm{d}\sigma(\boldsymbol{x}).

Invoking the trace inequality, using the fact $\|\,\nabla\operatorname{div}\boldsymbol{u}\,\|_{L^{2}}\leq\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\|_{L^{2}}$ and (2.10), we obtain

	$\displaystyle\iota^{2}\left\lvert\,\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{u}_{0}\mathrm{d}\sigma(\boldsymbol{x})\,\right\rvert$	$\displaystyle\leq C\iota^{2}\left(\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{H^{1}}^{1/2}\\|\,\boldsymbol{\epsilon}(\boldsymbol{u}_{0})\,\\|_{H^{1}}+\\|\,\nabla\operatorname{div}\boldsymbol{u}\,\\|_{L^{2}}^{1/2}\\|\,\nabla\operatorname{div}\boldsymbol{u}\,\\|_{H^{1}}^{1/2}\\|\,p_{0}\,\\|_{H^{1}}\right)$
		$\displaystyle\leq C\iota^{2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\boldsymbol{u}\,\\|_{H^{3}}^{1/2}\left(\\|\,\boldsymbol{u}_{0}\,\\|_{H^{2}}+\\|\,p_{0}\,\\|_{H^{1}}\right)$
		$\displaystyle\leq C\iota^{2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\boldsymbol{u}\,\\|_{H^{3}}^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}.$

Substituting (2.23) into (2.22), we obtain,

\|\,\boldsymbol{u}\,\|_{H^{3}}+\|\,p\,\|_{H^{2}}\leq C\iota^{-2}(\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}+\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}).

Invoking (2.26) again, we get

\iota^{2}\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\|_{L^{2}}^{1/2}\|\,\boldsymbol{u}\,\|_{H^{3}}^{1/2}\leq C\iota^{1/2}(\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}+\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}).

A combination of the above three inequalities gives

\iota^{2}\left\lvert\,\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{u}_{0}\mathrm{d}\sigma(\boldsymbol{x})\,\right\rvert\leq\dfrac{1}{4}\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}^{2}+C\iota\|\,\boldsymbol{f}\,\|_{L^{2}}^{2}.

Combining the above inequalities, we obtain

\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{\iota}\leq C\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}},

which together with (2.23) leads to

\|\,p-p_{0}\,\|_{\iota}\leq C\iota^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}.

The above two estimates immediately give (2.27) and

\|\,\nabla^{2}(\boldsymbol{u}-\boldsymbol{u}_{0})\,\|_{L^{2}}+\|\,\nabla(p-p_{0})\,\|_{L^{2}}\leq C\iota^{-1/2}\|\,\boldsymbol{f}\,\|_{L^{2}},

which together with (2.10) gives (2.28) with $k=0$ .

Substituting (2.27) into (2.22), we obtain the higher regularity estimate (2.28) with $k=1$ . ∎

A direct consequence of the above theorem is

Corollary 2.10.

There exists $C$ independent of $\iota$ and $\lambda$ such that

(2.29)

\|\,\boldsymbol{u}\,\|_{H^{3/2}}+\|\,p\,\|_{H^{1/2}}\leq C\|\,\boldsymbol{f}\,\|_{L^{2}},

and

(2.30)

\|\,\boldsymbol{u}\,\|_{H^{5/2}}+\|\,p\,\|_{H^{3/2}}\leq C\iota^{-1}\|\,\boldsymbol{f}\,\|_{L^{2}}.

Proof.

Using triangle inequality, (2.10) and (2.28), we obtain

\|\,\boldsymbol{u}-\boldsymbol{u}_{0}\,\|_{H^{2}}+\|\,p-p_{0}\,\|_{H^{1}}\leq\|\,\boldsymbol{u}\,\|_{H^{2}}+\|\,p\,\|_{H^{1}}+\|\,\boldsymbol{u}_{0}\,\|_{H^{2}}+\|\,p_{0}\,\|_{H^{1}}\leq C\iota^{-1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}.

Interpolating the above inequality and (2.27), we obtain

\|\,\boldsymbol{u}-\boldsymbol{u}_{0}\,\|_{H^{3/2}}+\|\,p-p_{0}\,\|_{H^{1/2}}\leq C\|\,\boldsymbol{f}\,\|_{L^{2}}.

Using (2.10), we get

\|\,\boldsymbol{u}_{0}\,\|_{H^{3/2}}+\|\,p_{0}\,\|_{H^{1/2}}\leq C\|\,\boldsymbol{u}_{0}\,\|_{H^{2}}+\|\,p_{0}\,\|_{H^{1}}\leq C\|\,\boldsymbol{f}\,\|_{L^{2}}.

A combination of the above two inequalities and the triangle inequality leads to (2.29).

Interpolating (2.28) with $k=0$ and (2.28) with $k=1$ , we obtain (2.30). ∎

3. A Family of Nonconforming Finite Elements

We introduce a family of finite elements to approximate the mixed variational problem (2.2). Let $\mathcal{T}_{h}$ be a triangulation of $\Omega$ with maximum mesh size $h$ . We assume all elements in $\mathcal{T}_{h}$ are shape-regular in the sense of Ciarlet and Raviart [Ciarlet:1978]. We also assume that $\mathcal{T}_{h}$ satisfies the inverse assumption: there exists $\sigma_{0}$ such that $h/h_{K}\leq\sigma_{0}$ for all $K\in\mathcal{T}_{h}$ . The space of piecewise vector fields is defined by

[H^{m}(\Omega,\mathcal{T}_{h})]^{2}{:}=\left\{\,\boldsymbol{v}\in[L^{2}(\Omega)]^{2}\,\mid\,\boldsymbol{v}|_{K}\in[H^{m}(T)]^{2}\quad\text{for all\quad}K\in\mathcal{T}_{h}\,\right\},

which is equipped with the broken norm

\|\,\boldsymbol{v}\,\|_{H^{m}(\Omega,\mathcal{T}_{h})}{:}=\|\,\boldsymbol{v}\,\|_{L^{2}}+\sum_{k=1}^{m}\|\,\nabla^{k}_{h}\boldsymbol{v}\,\|_{L^{2}},

where $\|\,\nabla^{k}_{h}\boldsymbol{v}\,\|_{L^{2}}^{2}=\sum_{K\in\mathcal{T}_{h}}\|\,\nabla^{k}\boldsymbol{v}\,\|_{L^{2}(K)}^{2}.$ For an interior edge $e$ shared by the triangles $K^{+}$ and $K^{-}$ , we define the jump of $\boldsymbol{v}$ across $e$ as

[\![\boldsymbol{v}]\!]{:}=\boldsymbol{v}^{+}\boldsymbol{n}^{+}+\boldsymbol{v}^{-}\boldsymbol{n}^{-}\qquad\text{with}\quad\boldsymbol{v}^{\pm}=v|_{K^{\pm}},

where $\boldsymbol{n}^{\pm}$ is the unit normal vector of $e$ towards the outside of $K^{\pm}$ . For $e\cap\partial\Omega\not=\emptyset$ , we set $[\![\boldsymbol{v}]\!]=\boldsymbol{v}\boldsymbol{n}$ .

3.1. A family of finite elements

Our construction is motivated by the element proposed in [GuzmanLeykekhmanNeilan:2012]. Define the element with a triple $(K,P_{K},\Sigma_{K})$ by specifying $K$ as a triangle, and

(3.1)

P_{K}{:}=\mathbb{P}_{r}(K)+b_{K}\sum_{i=1}^{3}b_{i}Q_{i}^{r-2}(K)+b_{K}^{2}R^{r-2}(K),

where $b_{K}=\prod_{i=1}^{3}\lambda_{i}$ and $b_{i}=b_{K}/\lambda_{i}$ with $\{\lambda_{i}\}_{i=1}^{3}$ the barycentric coordinates of $K$ .

Define

(3.2)

Q_{i}^{r-2}(K){:}=\left\{\,v\in\mathbb{P}_{r-2}(K)\,\mid\,\int_{K}b_{K}b_{i}vq\,\mathrm{d}\boldsymbol{x}=0\textup{ for all }q\in\mathbb{P}_{r-3}(K)\,\right\},

and

(3.3)

R^{r-2}(K){:}=\left\{\,v\in\mathbb{P}_{r-2}(K)\,\mid\,\int_{K}b_{K}^{2}vq\,\mathrm{d}\boldsymbol{x}=0\textup{ for all }q\in\mathbb{P}_{r-3}(K)\,\right\}.

The degrees of freedom (DoFs) for $P_{K}$ are given by

(3.4)

v\mapsto\left\{\begin{aligned} v(a)&\quad\text{for all vertices}\,a,\\ \int_{e}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{e}\partial_{\boldsymbol{n}}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{K}vq\,\mathrm{d}\boldsymbol{x}&\quad\text{for all\,}q\in\mathbb{P}_{r-2}(K).\end{aligned}\right.

We plot the DoFs for $r=2,3$ in Figure 1.

Figure 1. Diagram for DoFs. Left: DoFs for

r=2

are point evaluations of the function values at the vertex, the mean of the function along each edge, the mean of the normal derivative along each edge, and the mean of the function over the element; Right: DoFs for

r=3

are point evaluations of the function values at the vertex, the means of the function against

\mathbb{P}_{1}

along each edge, the means of the normal derivative against

\mathbb{P}_{1}

along each edge, and the means of the function against

\mathbb{P}_{1}

over the element

Lemma 3.1.

The set $(K,P_{K},\Sigma_{K})$ is unisolvent.

Proof.

Firstly we show that if the DoFs (3.4) can determine an element in $P_{K}$ , then the element is unique. Suppose $v\in P_{K}$ vanishes at the DoFs listed in (3.4), it suffices to show $v\equiv 0$ . Assume that

v=p_{r}+b_{K}\sum_{i=1}^{3}b_{i}q_{i}+b_{K}^{2}q_{r},

where $p_{r}\in\mathbb{P}_{r}(K)$ , and $q_{i}\in Q_{i}^{r-2}(K)$ , and $q_{r}\in R^{r-2}(K)$ . DoFs associated with $\mathbb{P}_{r}(K)$ are

v\mapsto\left\{\begin{aligned} v(a)&\quad\text{for all vertices}\,a,\\ \int_{e}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{K}vq\,\mathrm{d}\boldsymbol{x}&\quad\text{for all\,}q\in\mathbb{P}_{r-3}(K).\end{aligned}\right.

The bubble space vanishes on this subset of DoFs by (3.2) and (3.3). The number of the DoFs is $3+3(r-1)+(r-1)(r-2)/2=(r+1)(r+2)/2$ , which equals to the cardinality of $\mathbb{P}_{r}(K)$ . Hence $p_{r}\equiv 0$ .

A direct calculation gives

\int_{e_{i}}\partial_{\boldsymbol{n}}vq\,\mathrm{d}\sigma(\boldsymbol{x})=-\lvert\nabla\lambda_{i}\rvert\int_{e_{i}}b_{i}^{2}q_{i}q\,\mathrm{d}\sigma(\boldsymbol{x})=0\qquad\text{for all }q\in\mathbb{P}_{r-2}(e_{i}).

Taking $q=q_{i}$ in the above identity, we obtain $q_{i}=0$ on $e_{i}$ . Therefore, we write $q_{i}=\lambda_{i}p_{r-3}$ for certain $p_{r-3}\in\mathbb{P}_{r-3}(K)$ . Using (3.2), we get

\int_{K}b_{K}b_{i}q_{i}q\,\mathrm{d}\boldsymbol{x}=\int_{K}b_{K}^{2}p_{r-3}q\,\mathrm{d}\boldsymbol{x}=0\qquad\text{for all }q\in\mathbb{P}_{r-3}(K).

Taking $q=p_{r-3}$ in the above identity, we obtain $p_{r-3}\equiv 0$ . Therefore $v=b_{K}^{2}q_{r}$ for certain $q_{r}\in R^{r-2}(K)$ . The last set of DoFs equals zero, i.e.,

\int_{K}b_{K}^{2}q_{r}q\,\mathrm{d}\boldsymbol{x}=0\qquad\text{for all }q\in\mathbb{P}_{r-2}(K).

Taking $q=q_{r}$ in the above identity, we obtain $q_{r}\equiv 0$ . So does $v$ .

It remains to show the dimension of $P_{K}$ equals the number of DoFs (3.4). Proceeding along the same line as above, the element $v\equiv 0$ has a unique representation. Therefore (3.1) is a direct sum, and

\text{dim}P_{K}=\text{dim}\mathbb{P}_{r}(K)+4\left(\text{dim}\mathbb{P}_{r-2}(K)-\text{dim}\mathbb{P}_{r-3}(K)\right)=\dfrac{1}{2}(r^{2}+11r-6),

which equals to the number of DoFs (3.4) exactly. ∎

We define a local interpolation operator $\pi_{K}:H^{2}(K)\mapsto P_{K}$ as:

(3.5)

\left\{\begin{aligned} \pi_{K}v(a)&=v(a)&&\quad\text{for all vertices\,}a,\\ \int_{e}\pi_{K}vq\,\mathrm{d}\sigma(\boldsymbol{x})&=\int_{e}vq\,\mathrm{d}\sigma(\boldsymbol{x})&&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{e}\partial_{\boldsymbol{n}}\pi_{K}vq\,\mathrm{d}\sigma(\boldsymbol{x})&=\int_{e}\partial_{\boldsymbol{n}}vq\,\mathrm{d}\sigma(\boldsymbol{x})&&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{K}\pi_{K}vq\,\mathrm{d}\boldsymbol{x}&=\int_{K}vq\,\mathrm{d}\boldsymbol{x}&&\quad\text{for all\,}q\in\mathbb{P}_{r-2}(K).\end{aligned}\right.

Lemma 3.2.

There exists $C$ independent of $h_{K}$ such that for $v\in H^{k}(K)$ with $2\leq k\leq r+1$ , there holds

(3.6)

\|\,\nabla^{j}(v-\pi_{K}v)\,\|_{L^{2}(K)}\leq Ch_{K}^{k-j}\|\,\nabla^{k}v\,\|_{L^{2}(K)},

where $0\leq j\leq k$ .

Proof.

For any $v\in\mathbb{P}_{r}(K)\subset P_{K}$ , the definition (3.5) shows that $v-\pi_{K}v\in P_{K}$ and all DoFs of $v-P_{K}v$ vanish, then $v=\pi_{K}v$ . The estimate (3.6) immediately follows from the $\mathbb{P}_{r}(K)-$ invariance of the local interpolation operator $\pi_{K}$ [CiarletRaviart:1972]. ∎

Remark 3.3.

The element has a natural extension to three-dimensions by specifying $K$ as a tetrahedron, and

P_{K}{:}=\mathbb{P}_{r}(K)+b_{K}\sum_{i=1}^{4}b_{i}Q_{i}^{r-2}(K)+b_{K}^{2}R^{r-2}(K),

where $b_{K}=\prod_{i=1}^{4}\lambda_{i}$ is the element bubble function with $\lambda_{i}$ the barycentric coordinates associated with the vertices $a_{i}$ for $i=1,\cdots,4$ . $b_{i}=b_{K}/\lambda_{i}$ is the face bubble function associated with the face $f_{i}$ .

Define

Q_{i}^{r-2}(K){:}=\left\{\,v\in\mathbb{P}_{r-2}(K)\,\mid\,\int_{K}b_{K}b_{i}vq\,\mathrm{d}\boldsymbol{x}=0\textup{ for all }q\in\mathbb{P}_{r-3}(K)\,\right\},

and

R^{r-2}(K){:}=\left\{\,v\in\mathbb{P}_{r-2}(K)\,\mid\,\int_{K}b_{K}^{2}vq\,\mathrm{d}\boldsymbol{x}=0\textup{ for all }q\in\mathbb{P}_{r-4}(K)\,\right\}.

The DoFs for $P_{K}$ are given by

v\mapsto\left\{\begin{aligned} v(a)&\quad\text{for all vertices},\\ \int_{e}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all edges $e$ and\,}q\in\mathbb{P}_{r-2}(e),\\ \int_{f}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all faces $f$ and\,}q\in\mathbb{P}_{r-3}(f),\\ \int_{f}\partial_{\boldsymbol{n}}vq\,\mathrm{d}\sigma(\boldsymbol{x})&\quad\text{for all faces $f$ and\,}q\in\mathbb{P}_{r-2}(f),\\ \int_{K}vq\,\mathrm{d}\boldsymbol{x}&\quad\text{for all\,}q\in\mathbb{P}_{r-2}(K).\end{aligned}\right.

Similar to Lemma 3.1, the set $(K,P_{K},\Sigma_{K})$ is also unisolvent.

3.2. Explicit representation for the bubble space

We clarify the structures of (3.2) and (3.3) associated with the set of DoFs (3.4)₃ and the subset of (3.4)₄ respectively, and derive the explicit formulations of the corresponding shape functions, which seems missing in the literature, while such explicit representations are useful for implementation. We firstly recall the following facts about the Jacobi polynomials [Szego:1975]. For any $\alpha,\beta>-1$ and nonnegative integers $n,m$ , there holds

(3.7)

\int_{-1}^{1}(1-t)^{\alpha}(1+t)^{\beta}P_{n}^{(\alpha,\beta)}(t)P_{m}^{(\alpha,\beta)}(t)\mathrm{d}t=h_{n}^{(\alpha,\beta)}\delta_{nm},

where

h_{n}^{(\alpha,\beta)}=\int_{-1}^{1}(1-t)^{\alpha}(1+t)^{\beta}\bigl{[}P_{n}^{(\alpha,\beta)}(t)\bigr{]}^{2}\mathrm{d}t.

By [Szego:1975]*Eq. (4.3.3), we may write

(3.8)

h_{n}^{(\alpha,\beta)}=\dfrac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\dfrac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)\Gamma(n+1)},

where $\Gamma$ is the Gamma function.

One of the explicit form for $P_{n}^{(\alpha,\beta)}$ is

(1-t)^{\alpha}(1+t)^{\beta}P_{n}^{(\alpha,\beta)}(t)=\dfrac{(-1)^{n}}{2^{n}n!}\dfrac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\left((1-t)^{n+\alpha}(1+t)^{n+\beta}\right).

In particular,

P_{0}^{(\alpha,\beta)}(t)=1,\quad P_{1}^{(\alpha,\beta)}(t)=\dfrac{1}{2}(\alpha+\beta+2)t+\dfrac{1}{2}(\alpha-\beta).

Next we list certain facts about the Jacobi polynomials on the triangle [Dunkl:2014]*Section 2.4. For a triangle $K$ with vertices $a_{1},a_{2},a_{3}$ , any point $x\in K$ is uniquely expressed as

x=\lambda_{1}a_{1}+\lambda_{2}a_{2}+\lambda_{3}a_{3},\quad\lambda_{i}\geq 0\text{ and }\lambda_{1}+\lambda_{2}+\lambda_{3}=1.

Then $(\lambda_{1},\lambda_{2},\lambda_{3})$ is the barycentric coordinates of the point $x$ with respect to $K$ . For nonnegative integers $k,n$ such that $k\leq n$ , we define

(3.9)

P_{k,n}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3}){:}=(\lambda_{2}+\lambda_{3})^{k}P_{n-k}^{(2k+\beta+\gamma+1,\alpha)}(\lambda_{1}-\lambda_{2}-\lambda_{3})P_{k}^{(\gamma,\beta)}((\lambda_{2}-\lambda_{3})/(\lambda_{2}+\lambda_{3})).

It is straightforward to verify $P_{k,n}^{(\alpha,\beta,\gamma)}\in\mathbb{P}_{n}(K)$ . In particular,

P_{0,0}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})=1,\quad P_{0,1}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})=(\beta+\gamma+2)\lambda_{1}-(\alpha+1)(\lambda_{2}+\lambda_{3}),

and

P_{1,1}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})=(\gamma+1)\lambda_{2}-(\beta+1)\lambda_{3}.

For all $\alpha,\beta,\gamma>-1$ , and nonnegative integers $j,k,m,n$ such that $j\leq m$ and $k\leq n$ , there holds

(3.10)		$\displaystyle{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}\lambda_{1}^{\alpha}\lambda_{2}^{\beta}\lambda_{3}^{\gamma}P_{k,n}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})P_{j,m}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})\,\mathrm{d}\boldsymbol{x}$
	$\displaystyle=$	$\displaystyle 2\int_{\widehat{K}}\lambda_{1}^{\alpha}\lambda_{2}^{\beta}(1-\lambda_{1}-\lambda_{2})^{\gamma}P_{k,n}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},1-\lambda_{1}-\lambda_{2})P_{j,m}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},1-\lambda_{1}-\lambda_{2})\mathrm{d}\lambda_{1}\mathrm{d}\lambda_{2}$
	$\displaystyle=$	$\displaystyle 2h_{k,n}^{(\alpha,\beta,\gamma)}\delta_{jk}\delta_{mn},$

where $\widehat{K}:=\left\{\,(\lambda_{1},\lambda_{2})\,\mid\,\lambda_{1}\geq 0,\lambda_{2}\geq 0,\lambda_{1}+\lambda_{2}\leq 1\,\right\}$ is the standard reference triangle and

(3.11)	$\displaystyle h_{k,n}^{(\alpha,\beta,\gamma)}=$	$\displaystyle 2^{-(2k+\alpha+2\beta+2\gamma+3)}h_{n-k}^{(2k+\beta+\gamma+1,\alpha)}h_{k}^{(\gamma,\beta)}$
	$\displaystyle=$	$\displaystyle\dfrac{1}{(2n+\alpha+\beta+\gamma+2)(2k+\beta+\gamma+1)}$
		$\displaystyle\times\dfrac{\Gamma(n-k+\alpha+1)\Gamma(n+k+\beta+\gamma+2)\Gamma(k+\beta+1)\Gamma(k+\gamma+1)}{\Gamma(n-k+1)\Gamma(n+k+\alpha+\beta+\gamma+2)\Gamma(k+1)\Gamma(k+\beta+\gamma+1)}.$

Using the notation $(x)_{n}=\Gamma(x+n)/\Gamma(n)$ , we may find that the expression (3.11) is equivalent to the one in [Dunkl:2014]*Eq. (2.4.3). The identity (3.10) illustrates that $\{P_{k,n}^{(\alpha,\beta,\gamma)}(\lambda_{1},\lambda_{2},\lambda_{3})\mid 0\leq k\leq n\leq r\}$ are mutually orthogonal bases of $\mathbb{P}_{r}(K)$ with respect to the weight $\lambda_{1}^{\alpha}\lambda_{2}^{\beta}\lambda_{3}^{\gamma}$ .

Next we study the structure of the bubble spaces. For the barycentric coordinate function $\lambda_{i}$ such that $\lambda_{i}\equiv 0$ on $e_{i}$ , let $\lambda_{i}^{+}$ and $\lambda_{i}^{-}$ be the two other barycentric coordinates associated with the edges $e_{i}^{+}$ and $e_{i}^{-}$ , respectively. $(e_{i},e_{i}^{+},e_{i}^{-})$ are chosen in a counterclockwise manner. The space $Q_{i}^{r-2}(K)$ can be clarified by the Jacobi polynomials with respect to the weight $b_{K}b_{i}$ , while $R^{r-2}(K)$ can be clarified by the Jacobi polynomials with respect to the weight $b_{K}^{2}$ , which are formulated in the following lemmas.

Lemma 3.4.

The space $Q_{i}^{r-2}(K)$ takes the form

Q_{i}^{r-2}(K)=\textup{span}\left\{\,P_{k,r-2}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})\,\mid\,0\leq k\leq r-2\,\right\}.

Proof.

For any $v\in Q_{i}^{r-2}(K)\subset\mathbb{P}_{r-2}(K)$ , $v$ may be expanded into

v=\sum_{0\leq k\leq n\leq r-2}a_{kn}P_{k,n}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})

with unknown parameters $a_{kn}$ . Using the above representation, we may write the constraint in the definition (3.2) as

\sum_{0\leq k\leq n\leq r-2}a_{kn}{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}b_{K}b_{i}P_{k,n}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})q\mathrm{d}x=0\qquad\text{for all\quad}q\in\mathbb{P}_{r-3}(K).

Substituting $q=P_{j,m}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})$ for $0\leq j\leq m\leq r-3$ into the above equation, and using the orthogonal relation (3.10), we obtain $a_{jm}=0$ for $0\leq j\leq m\leq r-3$ . This concludes the lemma. ∎

Motivated by the above lemma, we change the definition of DoFs for the bubble space $b_{K}\sum_{i=1}^{3}b_{i}Q_{i}^{r-2}(K)$ from $\int_{e_{i}}\partial_{\boldsymbol{n}}vq\,\mathrm{d}\sigma(\boldsymbol{x})$ for any $q\in\mathbb{P}_{r-2}(e_{i})$ to

{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\dfrac{\partial v}{\partial n}P_{k}^{(2,2)}(\lambda_{i}^{+}-\lambda^{-})\,\mathrm{d}\sigma(\boldsymbol{x}),\quad k=0,\cdots,r-2.

Lemma 3.5.

The shape functions for the bubble space $b_{K}\sum_{i=1}^{3}b_{i}Q_{i}^{r-2}(K)$ associated with the above definition of DoFs are

a_{k,r-2}b_{K}b_{i}P_{k,r-2}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})

with

(3.12)

a_{k,r-2}=\dfrac{(-1)^{r-k-1}}{\lvert\nabla\lambda_{i}\rvert}\dfrac{(k+3)(k+4)(2k+5)}{(r-k-1)(k+1)(k+2)}\qquad k=0,\cdots,r-2.

Proof.

A direct calculation gives

\dfrac{\partial}{\partial\boldsymbol{n}}\left(b_{K}b_{i}P_{k,r-2}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})\right)|_{e_{i}}=-\lvert\nabla\lambda_{i}\rvert b_{i}^{2}P_{r-2-k}^{(2k+5,1)}(-1)P_{k}^{(2,2)}(\lambda_{i}^{+}-\lambda^{-}),

For $j=1,\cdots,r-2$ . Using the relation (3.7), and noting that $b_{i}=\lambda_{i}^{+}\lambda_{i}^{-}=\lambda_{i}^{+}(1-\lambda_{i}^{+})$ on $e_{i}$ , we obtain

		$\displaystyle{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\dfrac{\partial}{\partial\boldsymbol{n}}\left(b_{K}b_{i}P_{k,r-2}^{(1,2,2)}(\lambda_{i},\lambda_{i}^{+},\lambda_{i}^{-})\right)P_{j}^{(2,2)}(\lambda_{i}^{+}-\lambda^{-})\,\mathrm{d}\sigma(\boldsymbol{x})$
	$\displaystyle=$	$\displaystyle-\lvert\nabla\lambda_{i}\rvert P_{r-2-k}^{(2k+5,1)}(-1)\int_{0}^{1}\left(\lambda_{i}^{+}(1-\lambda_{i}^{+})\right)^{2}P_{k}^{(2,2)}(2\lambda_{i}^{+}-1)P_{j}^{(2,2)}(2\lambda_{i}^{+}-1)\mathrm{d}\lambda_{i}^{+}$
	$\displaystyle=$	$\displaystyle-\dfrac{\lvert\nabla\lambda_{i}\rvert}{32}P_{r-2-k}^{(2k+5,1)}(-1)\int_{-1}^{1}(1-t)^{2}(1+t)^{2}P_{k}^{(2,2)}(t)P_{j}^{(2,2)}(t)\mathrm{d}t$
	$\displaystyle=$	$\displaystyle-\dfrac{\lvert\nabla\lambda_{i}\rvert}{32}P_{r-2-k}^{(2k+5,1)}(-1)h_{k}^{(2,2)}\delta_{jk}.$

This gives

a_{k,r-2}=-\dfrac{32}{\lvert\nabla\lambda_{i}\rvert P_{r-2-k}^{(2k+5,1)}(-1)h_{k}^{(2,2)}}.

By [Szego:1975]*Eq. (4.1.1),(4.1.3), we obtain $P_{r-2-k}^{(2k+5,1)}(-1)=(-1)^{r-k}(r-k-1).$ Using (3.8), we obtain

h_{k}^{(2,2)}=\dfrac{32(k+1)(k+2)}{(k+3)(k+4)(2k+5)}.

A combination of the above three identities leads to (3.12). ∎

Next we list the shape functions for the elements of low-order.

Example 3.6.

The bubble space $b_{K}\sum_{i=1}^{3}b_{i}Q_{i}^{r-2}(K)$ for the lowest-order $r=2$ is

b_{K}\textup{span}\left\{\,b_{i}\,\mid\,i=1,2,3\,\right\}.

The shape functions associated with ${\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\partial_{\boldsymbol{n}}v\,\mathrm{d}\sigma(\boldsymbol{x})$ is $-\dfrac{30}{\lvert\nabla\lambda_{i}\rvert}b_{K}b_{i}.$

The bubble space $b_{K}\sum_{i=1}^{3}b_{i}Q_{i}^{r-2}(K)$ for the case $r=3$ is

b_{K}\textup{span}\left\{\,b_{i}(3\lambda_{i}-\lambda_{i}^{+}-\lambda_{i}^{-}),b_{i}(\lambda_{i}^{+}-\lambda_{i}^{-})\,\mid\,i=1,2,3\,\right\}.

The shape functions associated with ${\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\partial_{\boldsymbol{n}}v\,\mathrm{d}\sigma(\boldsymbol{x})$ is

\dfrac{30}{\lvert\nabla\lambda_{i}\rvert}b_{K}b_{i}(3\lambda_{i}-\lambda_{i}^{+}-\lambda_{i}^{-}).

The shape functions associated with $3{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\partial_{\boldsymbol{n}}v(\lambda_{i}^{+}-\lambda_{i}^{-})\,\mathrm{d}\sigma(\boldsymbol{x})$ is

-\dfrac{70}{\lvert\nabla\lambda_{i}\rvert}b_{K}b_{i}(\lambda_{i}^{+}-\lambda_{i}^{-}).

Lemma 3.7.

The space $R^{r-2}(K)$ takes the forms

R^{r-2}(K)=\textup{span}\left\{\,P_{k,r-2}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})\,\mid\,0\leq k\leq r-2\,\right\}.

Proof.

For any $v\in R^{r-2}(K)\subset\mathbb{P}_{r-2}(K)$ , we expand $v$ into

v=\sum_{0\leq k\leq n\leq r-2}a_{kn}P_{k,n}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})

with $a_{kn}$ to be determined later on. Using the above representation, we may write the constraint in the definition of the space (3.3) as: For all $q\in\mathbb{P}_{r-3}(K)$ ,

\sum_{0\leq k\leq n\leq r-2}a_{kn}{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}b_{K}^{2}P_{k,n}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})q\mathrm{d}x=0.

We substitute $q=P_{j,m}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})$ for $0\leq j\leq m\leq r-3$ into the above equation. Using the orthogonal relation (3.10), we may obtain $a_{jm}=0$ for $0\leq j\leq m\leq r-3$ . This completes the proof. ∎

Motivated by the above lemma, we change the definition of DoFs for $b_{K}^{2}R^{r-2}(K)$ from $\int_{K}vq\,\mathrm{d}\boldsymbol{x}$ for any $q\in\mathbb{P}_{r-2}(K)\setminus\mathbb{P}_{r-3}(K)$ to

{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}vP_{k,r-2}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})\,\mathrm{d}\boldsymbol{x},\quad k=0,\cdots,r-2.

Lemma 3.8.

The shape functions for the bubble space $b_{K}^{2}R^{r-2}(K)$ associated with the above definition of DoFs are

a_{k,r-2}b_{K}^{2}P_{k,r-2}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})

with

a_{k,r-2}=\dfrac{(2r+4)(2k+5)(r+k+4)(r+k+5)(k+3)(k+4)}{2(r-k)(r-1-k)(k+1)(k+2)}\qquad k=0,\cdots,r-2.

Proof.

For $j=1,\cdots,r-2$ , we obtain

{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}b_{K}^{2}P_{k,r-2}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})P_{j,r-2}^{(2,2,2)}(\lambda_{1},\lambda_{2},\lambda_{3})\,\mathrm{d}\boldsymbol{x}=2h_{k,r-2}^{(2,2,2)}\delta_{jk},

which gives

a_{k,r-2}=\dfrac{1}{2h_{k,r-2}^{(2,2,2)}}.

Using (3.11), we obtain

h_{k,r-2}^{(2,2,2)}=\dfrac{(r-k)(r-1-k)(k+1)(k+2)}{(2r+4)(2k+5)(r+k+4)(r+k+5)(k+3)(k+4)}.

These give the simplified expression of $a_{k,r-2}$ . ∎

According to the definition of bubble space, we may have

Example 3.9.

The bubble space $b_{K}^{2}R^{r-2}(K)$ for the lowest-order case $r=2$ is span $\{b_{K}^{2}\}$ .

The shape functions associated with ${\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}v\,\mathrm{d}\boldsymbol{x}$ is $2520b_{K}^{2}$ .

The bubble space $b_{K}^{2}R^{r-2}(K)$ for the case $r=3$ is $b_{K}^{2}$ span $\{2\lambda_{1}-\lambda_{2}-\lambda_{3},\lambda_{2}-\lambda_{3}\}$ .

The shape functions associated with $3\,{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}v(2\lambda_{1}-\lambda_{2}-\lambda_{3})\,\mathrm{d}\boldsymbol{x}$ is $4200b_{K}^{2}(2\lambda_{1}-\lambda_{2}-\lambda_{3})$ .

The shape functions associated with $3\,{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}v(\lambda_{2}-\lambda_{3})\,\mathrm{d}\boldsymbol{x}$ is $12600b_{K}^{2}(\lambda_{2}-\lambda_{3})$ .

Remark 3.10.

Based on the above results, we give the details for the element of the lowest-order, i.e., $r=2$ , which have been used in the numerical examples. The local finite element space

P_{K}=\mathbb{P}_{2}(K)+b_{K}\sum_{i=1}^{3}\text{span}\{b_{i}\}+\text{span}\{b_{K}^{2}\},

and DoFs

\Sigma_{K}=\left\{\,v(a_{i}),{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}v\,\mathrm{d}\sigma(\boldsymbol{x}),{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\partial_{\boldsymbol{n}}v\,\mathrm{d}\sigma(\boldsymbol{x}),{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}v\,\mathrm{d}\boldsymbol{x}\,\mid\,i=1,2,3\,\right\}.

The shape functions associated with $\{v(a_{i})\}_{i=1,2,3}$ are

\phi_{i}=\lambda_{i}(3\lambda_{i}-2)+30b_{K}\left(2b_{i}+\sum_{j\neq i}\dfrac{\nabla\lambda_{i}\cdot\nabla\lambda_{j}}{\lvert\nabla\lambda_{j}\rvert^{2}}b_{j}(4\lambda_{j}-1)+6b_{K}\right).

The shape functions associated with $\{{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}v\mathrm{d}\sigma(\boldsymbol{x})\}_{i=1,2,3}$ are

\varphi_{i}=6b_{i}+90b_{K}\left(b_{i}-\sum_{j\neq i}b_{j}-10b_{K}\right).

The shape functions associated with $\{{\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{e_{i}}\partial_{\boldsymbol{n}}v\mathrm{d}\sigma(\boldsymbol{x})\}_{i=1,2,3}$ are

\psi_{i}=\dfrac{30}{\lvert\nabla\lambda_{i}\rvert}b_{K}b_{i}(4\lambda_{i}-1).

The shape functions associated with ${\int\negthickspace\negthickspace\negthickspace\negthinspace-}_{K}v\,\mathrm{d}\boldsymbol{x}$ is $\phi_{0}=2520b_{K}^{2}$ .

4. The Mixed Finite Elements Approximation

In this part, we construct stable finite element pairs to approximate (2.2). We ignore the influence of the curved boundary for error estimates for brevity. Define

X_{h}{:}=\left\{\,v\in H_{0}^{1}(\Omega)\,\mid\,v|_{K}\in P_{K}\text{ for all }K\in\mathcal{T}_{h},\quad\int_{e}[\![\partial_{\boldsymbol{n}}v]\!]q\,\mathrm{d}\sigma(\boldsymbol{x})=0\text{ for all }e\in\mathcal{E}_{h},q\in\mathbb{P}_{r-2}(e)\,\right\},

and $V_{h}{:}=[X_{h}]^{2}$ . Let $P_{h}\subset P$ be the continuous Lagrangian finite element of order $r-1$ . We shall prove a uniform discrete B-B inequality for the pair $(V_{h},P_{h})$ .

The following rescaled trace inequality will be used later on: There exists $C$ independent of $h_{K}$ such that

(4.1)

\|\,v\,\|_{L^{2}(\partial K)}\leq C\left(h_{K}^{-1/2}\|\,v\,\|_{L^{2}(K)}+\|\,v\,\|_{L^{2}(K)}^{1/2}\|\,\nabla\boldsymbol{v}\,\|_{L^{2}(K)}^{1/2}\right).

This inequality may be found in [BrennerScott:2008].

4.1. The mixed finite element approximation

We define the mixed finite element approximation problem as follows. Find $\boldsymbol{u}_{h}\in V_{h}$ and $p_{h}\in P_{h}$ such that

(4.2)

\left\{\begin{aligned} a_{\iota,h}(\boldsymbol{u}_{h},\boldsymbol{v})+b_{\iota,h}(\boldsymbol{v},p_{h})&=(\boldsymbol{f},\boldsymbol{v})\qquad&&\text{for all\quad}\boldsymbol{v}\in V_{h},\\ b_{\iota,h}(\boldsymbol{u}_{h},q)-\lambda^{-1}c_{\iota}(p_{h},q)&=0\qquad&&\text{for all\quad}q\in P_{h},\end{aligned}\right.

where

	$\displaystyle a_{\iota,h}(\boldsymbol{v},\boldsymbol{w}){:}$	$\displaystyle=2\mu\Bigl{(}(\boldsymbol{\epsilon}(\boldsymbol{v}),\boldsymbol{\epsilon}(\boldsymbol{w}))+\iota^{2}(\nabla_{h}\boldsymbol{\epsilon}(\boldsymbol{v}),\nabla_{h}\boldsymbol{\epsilon}(\boldsymbol{w}))\Bigr{)}\quad$		$\displaystyle\text{for all\,}\boldsymbol{v},\boldsymbol{w}\in V_{h},$
	$\displaystyle b_{\iota,h}(\boldsymbol{v},q){:}$	$\displaystyle=(\operatorname{div}\boldsymbol{v},q)+\iota^{2}(\nabla_{h}\operatorname{div}\boldsymbol{v},\nabla q)\quad$		$\displaystyle\text{for all\,}\boldsymbol{v}\in V_{h},q\in P_{h}.$

Note that $V_{h}\not\subset V$ , and this is a nonconforming method, we introduce the broken norm

\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}{:}=\|\,\nabla\boldsymbol{v}\,\|_{L^{2}}+\iota\|\,\nabla_{h}^{2}\boldsymbol{v}\,\|_{L^{2}}\qquad\text{for all}\quad\boldsymbol{v}\in V_{h}

Due to the continuity of $\boldsymbol{v}$ , $\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}$ is a norm over $V_{h}$ .

The following broken Korn’s inequality was proved in [LiMingWang:2021]*Theorem 2:

\|\,\nabla_{h}\boldsymbol{\epsilon}(\boldsymbol{v})\,\|_{L^{2}}\geq\left(1-1/\sqrt{2}\right)\|\,\nabla_{h}^{2}\boldsymbol{v}\,\|_{L^{2}},

which together with the first Korn’s inequality (2.6) gives

(4.3)

a_{\iota,h}(\boldsymbol{v},\boldsymbol{v})\geq\dfrac{\mu}{2}\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}^{2}\qquad\text{for all }\boldsymbol{v}\in V_{h}.

It remains to prove the discrete B-B inequality for the pair $(V_{h},P_{h})$ . To this end, we construct a Fortin operator that is uniformly stable in the weighted norm $\|\,\nabla\cdot\,\|_{\iota,h}$ [Winther:2013]. The key is to construct different Fortin operators for $\iota/h$ in different regimes.

Firstly we define an interpolation operator $\varPi_{h}:V\to V_{h}$ by $\varPi_{h}|_{K}{:}=\varPi_{K}=[\pi_{K}]^{2}$ , which satisfies

Lemma 4.1.

For all $\boldsymbol{v}\in V$ , there holds

(4.4)

b_{\iota,h}(\varPi_{h}\boldsymbol{v},p)=b_{\iota}(\boldsymbol{v},p)\qquad\text{for all\;}p\in P_{h}.

Proof.

Using the fact that $\varPi_{h}\boldsymbol{v}\in V_{h}\subset[H_{0}^{1}(\Omega)]^{2}$ , an integration by parts gives

(4.5)

\int_{\Omega}\operatorname{div}(\boldsymbol{v}-\varPi_{h}\boldsymbol{v})p\,\mathrm{d}\boldsymbol{x}=-\sum_{K\in\mathcal{T}_{h}}\int_{K}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\cdot\nabla p\,\mathrm{d}\boldsymbol{x}=0,

where we have used the identity (3.5)₄ in the last step.

Next, integration by parts yields

\int_{\Omega}\nabla\operatorname{div}(\boldsymbol{v}-\varPi_{h}\boldsymbol{v})\cdot\nabla p\,\mathrm{d}\boldsymbol{x}=\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}\operatorname{div}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\partial_{\boldsymbol{n}}p\,\mathrm{d}\sigma(\boldsymbol{x})-\sum_{K\in\mathcal{T}_{h}}\int_{K}\operatorname{div}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\Delta p\,\mathrm{d}\boldsymbol{x}=0.

The first term vanishes because $\partial_{j}=t_{j}\partial_{\boldsymbol{t}}+n_{j}\partial_{\boldsymbol{n}}$ for components $j=1,2$ , and using (3.5)₂, we obtain that for each edge $e\in\partial K$ ,

\int_{e}t_{j}\dfrac{\partial}{\partial\boldsymbol{t}}(v_{j}-\pi_{K}v_{j})\partial_{\boldsymbol{n}}p\,\mathrm{d}\sigma(\boldsymbol{x})=-\int_{e}t_{j}(v_{j}-\pi_{K}v_{j})\dfrac{\partial^{2}p}{\partial t\partial n}\,\mathrm{d}\sigma(\boldsymbol{x})=0.

Using (3.5)₃, we obtain

\int_{e}n_{j}\dfrac{\partial}{\partial\boldsymbol{n}}(v_{j}-\pi_{K}v_{j})\partial_{\boldsymbol{n}}p\,\mathrm{d}\sigma(\boldsymbol{x})=0.

While the second term vanishes because

\int_{K}\operatorname{div}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\Delta p\,\mathrm{d}\boldsymbol{x}=\int_{\partial K}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\cdot\boldsymbol{n}\Delta p\,\mathrm{d}\sigma(\boldsymbol{x})-\int_{K}(\boldsymbol{v}-\varPi_{K}\boldsymbol{v})\cdot\nabla\Delta p\,\mathrm{d}\boldsymbol{x}=0,

where we have used (3.5)₂ and (3.5)₄. ∎

The operator $\varPi_{h}$ is not H ${}^{1}-$ bounded by (3.6), and we construct an H ${}^{1}-$ bounded regularized interpolation operator as follows.

Lemma 4.2.

There exists an operator $I_{h}:\,V\mapsto V_{h}$ satisfying

(4.6)

\int_{\Omega}\operatorname{div}(\boldsymbol{v}-I_{h}\boldsymbol{v})p\,\mathrm{d}\boldsymbol{x}=0\qquad\text{for all\quad}p\in P_{h},

and if $\boldsymbol{v}\in V\cap[H^{s}(\Omega)]^{2}$ with $1\leq s\leq r+1$ , then

(4.7)

\|\,\nabla^{j}_{h}(\boldsymbol{v}-I_{h}\boldsymbol{v})\,\|_{L^{2}}\leq Ch^{s-j}\|\,\nabla^{s}\boldsymbol{v}\,\|_{L^{2}}\qquad 0\leq j\leq s.

The construction of $I_{h}$ is based on a regularized interpolation operator in [GuzmanLeykekhmanNeilan:2012] and the standard construction of the Fortin operator [Fortin:1977]. The operator $I_{h}$ is also well-defined for functions in $[H^{2}(\Omega)\cap H_{0}^{1}(\Omega)]^{2}$ .

Proof.

Define $I_{h}\,:V\mapsto V_{h}$ with $I_{h}{:}=[\varPi_{1}]^{2}$ and

\varPi_{1}{:}=\varPi_{0}(I-\varPi_{2})+\varPi_{2},

where the regularized interpolation operator $\varPi_{2}:H_{0}^{2}(\Omega)\mapsto X_{h}$ was constructed in [GuzmanLeykekhmanNeilan:2012]*Lemma 2, which satisfies

(4.8)

\|\,\nabla^{j}_{h}(v-\varPi_{2}v)\,\|_{L^{2}}\leq Ch^{s-j}\|\,\nabla^{s}v\,\|_{L^{2}},\qquad 1\leq s\leq r+1,0\leq j\leq s.

The operator $\varPi_{0}\,:H_{0}^{1}(\Omega)\mapsto X_{h}$ is defined for any element $K\in\mathcal{T}_{h}$ as

\left\{\begin{aligned} \varPi_{0}v(a)&=0\quad\text{for all vertices}\,a,\\ \int_{e}\varPi_{0}vq\,\mathrm{d}\sigma(\boldsymbol{x})&=0\quad\forall q\in\mathbb{P}_{r-2}(e)\quad\text{for all edges}\,e,\\ \int_{e}\partial_{\boldsymbol{n}}\varPi_{0}vq\,\mathrm{d}\sigma(\boldsymbol{x})&=0\quad\forall q\in\mathbb{P}_{r-2}(e)\quad\text{for all edges}\,e,\\ \int_{K}\varPi_{0}vq\,\mathrm{d}\boldsymbol{x}&=\int_{K}vq\,\mathrm{d}\boldsymbol{x}\quad\text{for all\quad}q\in\mathbb{P}_{r-2}(K).\end{aligned}\right.

On each element $K$ , we have $(\varPi_{0}v)|_{K}\in H_{0}^{1}(K)$ for any $v\in H_{0}^{1}(\Omega)$ , and hence $\varPi_{0}v\in X_{h}$ . A standard scaling argument gives

\|\,\varPi_{0}v\,\|_{L^{2}}\leq C\|\,v\,\|_{L^{2}}\qquad\text{for all }v\in H_{0}^{1}(\Omega).

For any $\boldsymbol{v}\in V$ , using the fact that $I_{h}\boldsymbol{v}\in V_{h}$ and $p\in P_{h}$ , an integration by parts gives

	$\displaystyle\int_{\Omega}\operatorname{div}(\boldsymbol{v}-I_{h}\boldsymbol{v})p\,\mathrm{d}\boldsymbol{x}$	$\displaystyle=-\sum_{K\in\mathcal{T}_{h}}\int_{K}(\boldsymbol{v}-I_{h}\boldsymbol{v})\cdot\nabla p\,\mathrm{d}\boldsymbol{x}+\sum_{K\in\mathcal{T}_{h}}\int_{\partial K}(\boldsymbol{v}-I_{h}\boldsymbol{v})\cdot\boldsymbol{n}p\,\mathrm{d}\sigma(\boldsymbol{x})$
		$\displaystyle=-\sum_{K\in\mathcal{T}_{h}}\int_{K}(\boldsymbol{v}-I_{h}\boldsymbol{v})\cdot\nabla p\,\mathrm{d}\boldsymbol{x}=-\sum_{K\in\mathcal{T}_{h}}\int_{K}(I-\varPi_{0})(I-\varPi_{2})v\cdot\nabla p\,\mathrm{d}\boldsymbol{x}$
		$\displaystyle=0,$

where we have used the last property of $\varPi_{0}$ . This gives (4.6).

Using (4.8), the L²-stability of $\varPi_{0}$ and the inverse inequality, we obtain

	$\displaystyle\\|\,\nabla^{j}_{h}(\boldsymbol{v}-\varPi_{1}\boldsymbol{v})\,\\|_{L^{2}}\leq$	$\displaystyle\\|\,\nabla^{j}_{h}(\boldsymbol{v}-\varPi_{2}\boldsymbol{v})\,\\|_{L^{2}}+\\|\,\nabla^{j}_{h}\varPi_{0}(\boldsymbol{v}-\varPi_{2}\boldsymbol{v})\,\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{s-j}\\|\,\nabla^{s}\boldsymbol{v}\,\\|_{L^{2}}+Ch^{-j}\\|\,\boldsymbol{v}-\varPi_{2}\boldsymbol{v}\,\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{s-j}\\|\,\nabla^{s}\boldsymbol{v}\,\\|_{L^{2}}.$

This implies (4.7) and completes the proof. ∎

We are ready to prove the following discrete B-B inequality.

Theorem 4.3.

There exists $\beta$ independent of $\iota$ and $h$ , such that

(4.9)

\sup_{\boldsymbol{v}\in V_{h}}\dfrac{b_{\iota,h}(\boldsymbol{v},p)}{\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}}\geq\beta\|p\|_{\iota}\qquad\text{for all\quad}p\in P_{h}.

Proof.

Using (2.3), for any $p\in P_{h}\subset P$ , there exists $\boldsymbol{v}_{0}\in V$ such that

b_{\iota}(\boldsymbol{v}_{0},p)=\|p\|_{\iota}^{2}\qquad\text{and\quad}\|\nabla\boldsymbol{v}_{0}\|_{\iota}\leq C\|p\|_{\iota}.

First, we consider the case $\iota/h\leq\gamma$ with $\gamma$ to be determined later on. By (4.7), we obtain

\|\,\nabla I_{h}\boldsymbol{v}_{0}\,\|_{\iota,h}\leq\|\nabla\boldsymbol{v}_{0}\|_{\iota}+\|\,\nabla(\boldsymbol{v}_{0}-I_{h}\boldsymbol{v}_{0})\,\|_{\iota,h}\leq C\|\nabla\boldsymbol{v}_{0}\|_{\iota},

and

\|\,\nabla_{h}\operatorname{div}(\boldsymbol{v}_{0}-I_{h}\boldsymbol{v}_{0})\,\|_{L^{2}}\leq C\|\,\nabla^{2}\boldsymbol{v}_{0}\,\|_{L^{2}}\leq C\|\,\nabla p\,\|_{L^{2}}.

Combining the above inequality and using the inverse inequality for any $p\in P_{h}$ , we obtain

\iota^{2}\lvert(\nabla_{h}\operatorname{div}(\boldsymbol{v}_{0}-I_{h}\boldsymbol{v}_{0}),\nabla p)\rvert\leq C\iota^{2}\|\,\nabla p\,\|_{L^{2}}^{2}\leq C_{\ast}(\iota/h)^{2}\|\,p\,\|_{L^{2}}^{2}\leq\gamma^{2}C_{\ast}\|p\|_{\iota}^{2}.

Fix $\gamma$ such that $\gamma^{2}C_{\ast}<1$ , we obtain

b_{\iota,h}(I_{h}\boldsymbol{v}_{0},p)=b_{\iota}(\boldsymbol{v}_{0},p)-\iota^{2}(\nabla_{h}\operatorname{div}(\boldsymbol{v}_{0}-I_{h}\boldsymbol{v}_{0}),\nabla p)\geq(1-\gamma^{2}C_{\ast})\|p\|_{\iota}^{2}.

This gives

(4.10)

\sup_{v\in V_{h}}\dfrac{b_{\iota,h}(\boldsymbol{v},p)}{\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}}\geq\dfrac{b_{\iota,h}(I_{h}\boldsymbol{v}_{0},p)}{\|\,\nabla I_{h}v_{0}\,\|_{\iota,h}}\geq\dfrac{1-\gamma^{2}C_{*}}{C}\|p\|_{\iota}.

Next, if $\iota/h>\gamma$ , then we use (3.6) and obtain

\|\,\nabla(\boldsymbol{v}_{0}-\varPi_{h}\boldsymbol{v}_{0})\,\|_{L^{2}}\leq Ch\|\,\nabla^{2}\boldsymbol{v}_{0}\,\|_{L^{2}},\quad\text{and}\quad\|\,\nabla^{2}_{h}(\boldsymbol{v}_{0}-\varPi_{h}\boldsymbol{v}_{0})\,\|_{L^{2}}\leq C\|\,\nabla^{2}\boldsymbol{v}_{0}\,\|_{L^{2}}.

Therefore,

\|\,\nabla(\boldsymbol{v}_{0}-\varPi_{h}\boldsymbol{v}_{0})\,\|_{\iota,h}\leq C(h+\iota)\|\,\nabla^{2}\boldsymbol{v}_{0}\,\|_{L^{2}}\leq C(1+h/\iota)\|\nabla\boldsymbol{v}_{0}\|_{\iota}\leq C(1+1/\gamma)\|\nabla\boldsymbol{v}_{0}\|_{\iota}.

Hence,

\|\,\nabla\varPi_{h}\boldsymbol{v}_{0}\,\|_{\iota,h}\leq\|\nabla\boldsymbol{v}_{0}\|_{\iota}+\|\,\nabla(\boldsymbol{v}_{0}-\varPi_{h}\boldsymbol{v}_{0})\,\|_{\iota,h}\leq C(2+1/\gamma)\|\nabla\boldsymbol{v}_{0}\|_{\iota},

which together with (4.4) gives

(4.11)

\sup_{v\in V_{h}}\dfrac{b_{\iota,h}(\boldsymbol{v},p)}{\|\,\nabla\boldsymbol{v}\,\|_{\iota,h}}\geq\dfrac{b_{\iota,h}(\varPi_{h}\boldsymbol{v}_{0},p)}{\|\,\nabla\varPi_{h}\boldsymbol{v}_{0}\,\|_{\iota,h}}=\dfrac{b_{\iota,h}(\boldsymbol{v}_{0},p)}{\|\,\nabla\varPi_{h}\boldsymbol{v}_{0}\,\|_{\iota,h}}=\dfrac{\|p\|_{\iota}^{2}}{\|\,\nabla\varPi_{h}\boldsymbol{v}_{0}\,\|_{\iota,h}}\geq\dfrac{\gamma}{C(1+2\gamma)}\|p\|_{\iota}.

A combination of (4.10) and (4.11) shows that (4.9) holds true with $\beta$ independent of $\iota$ and $h$ . ∎

The well-posedness of the mixed approximation problem (4.2) follows from the ellipticity of $a_{\iota,h}$ and the discrete B-B inequality of $b_{\iota,h}$ . We are ready to derive the error estimate.

4.2. Error estimates

To carry out the error estimate, we define the bilinear form $\mathcal{A}$ as

\mathcal{A}(\boldsymbol{v},q;\boldsymbol{w},z){:}=a_{\iota,h}(\boldsymbol{v},\boldsymbol{w})+b_{\iota,h}(\boldsymbol{w},q)+b_{\iota,h}(\boldsymbol{v},z)-\lambda^{-1}c_{\iota}(q,z)

for all $\boldsymbol{v},\boldsymbol{w}\in V_{h}$ and $q,z\in P_{h}$ .

We prove the following inf-sup inequality for $\mathcal{A}$ with the aid of the discrete B-B inequality (4.9).

Lemma 4.4.

There exists $\alpha$ depending on $\mu$ and $\beta$ such that

(4.12)

\inf_{(\boldsymbol{v},q)\in V_{h}\times P_{h}}\sup_{(\boldsymbol{w},z)\in V_{h}\times P_{h}}\dfrac{\mathcal{A}(\boldsymbol{v},q;\boldsymbol{w},z)}{|\!|\!|(\boldsymbol{w},z)|\!|\!|_{\iota,h}|\!|\!|(\boldsymbol{v},q)|\!|\!|_{\iota,h}}\geq\alpha,

where $|\!|\!|(\boldsymbol{w},z)|\!|\!|_{\iota,h}{:}=\|\,\nabla\boldsymbol{w}\,\|_{\iota,h}+\|z\|_{\iota}+\lambda^{-1/2}\|z\|_{\iota}$ and $\beta$ has appeared in (4.9).

Proof.

Noting that $a_{\iota,h}$ is elliptic over $V_{h}$ (4.3) and the discrete B-B inequality for $b_{\iota,h}$ holds (4.9), we obtain (4.12) by [Braess:1996]*Theorem 2. ∎

We are ready to prove error estimates.

Theorem 4.5.

There exists $C$ independent of $\iota,\lambda$ and $h$ such that

(4.13)

|\!|\!|(\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h})|\!|\!|_{\iota,h}\leq C(h^{r}+\iota h^{r-1})(\|\,\boldsymbol{u}\,\|_{H^{r+1}}+\|\,p\,\|_{H^{r}}),

and

(4.14)

|\!|\!|(\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h})|\!|\!|_{\iota,h}\leq Ch^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}}.

Proof.

Let $\boldsymbol{v}=\boldsymbol{u}_{h}-\boldsymbol{u}_{I}$ and $q=p_{h}-p_{I}$ with $\boldsymbol{u}_{I}\in V_{h}$ and $p_{I}\in P_{h}$ , for any $\boldsymbol{w}\in V_{h}$ and $z\in P_{h}$ ,

	$\displaystyle\mathcal{A}(\boldsymbol{v},q;\boldsymbol{w},z)$	$\displaystyle=\mathcal{A}(\boldsymbol{u}_{h},p_{h};\boldsymbol{w},z)-\mathcal{A}(\boldsymbol{u},p;\boldsymbol{w},z)+\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I};\boldsymbol{w},z)$
		$\displaystyle=(\boldsymbol{f},\boldsymbol{w})-\mathcal{A}(\boldsymbol{u},p;\boldsymbol{w},z)+\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I};\boldsymbol{w},z)$
		$\displaystyle=\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I};\boldsymbol{w},z)-\iota^{2}\sum_{e\in\mathcal{E}_{h}}\int_{e}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot[\![\partial_{\boldsymbol{n}}\boldsymbol{w}]\!]\,\mathrm{d}\sigma(\boldsymbol{x}).$

The boundedness of $\mathcal{A}$ yields

\lvert\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I};\boldsymbol{w},z)\rvert\leq\max(1,2\mu)|\!|\!|(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I})|\!|\!|_{\iota,h}|\!|\!|(\boldsymbol{w},z)|\!|\!|_{\iota,h}.

Let $\boldsymbol{u}_{I}=\varPi_{h}\boldsymbol{u}$ be the interpolation of $\boldsymbol{u}$ and $p_{I}$ be the $r-1$ order Lagrangian interpolation of $p$ , respectively. The standard interpolation error estimates in (3.6) gives

|\!|\!|(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I})|\!|\!|_{\iota,h}\leq C(h^{r}+\iota h^{r-1})(\|\,\boldsymbol{u}\,\|_{H^{r+1}}+\|\,p\,\|_{H^{r}}).

Note that

\int_{e}[\![\partial_{\boldsymbol{n}}\boldsymbol{w}]\!]q\,\mathrm{d}\sigma(\boldsymbol{x})=0\qquad\text{for all }q\in\mathbb{P}_{r-2}(e).

A standard estimate for the consistency error functional with trace inequality (4.1) gives

	$\displaystyle\iota^{2}\left\lvert\,\sum_{e\in\mathcal{E}_{h}}\int_{e}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot[\![\partial_{\boldsymbol{n}}\boldsymbol{w}]\!]\,\mathrm{d}\sigma(\boldsymbol{x})\,\right\rvert$	$\displaystyle\leq C\iota^{2}h^{r-1}(\\|\,\boldsymbol{u}\,\\|_{H^{r+1}}+\\|\,p\,\\|_{H^{r}})\\|\,\nabla_{h}^{2}\boldsymbol{w}\,\\|_{L^{2}}$
		$\displaystyle\leq C\iota h^{r-1}(\\|\,\boldsymbol{u}\,\\|_{H^{r+1}}+\\|\,p\,\\|_{H^{r}})\\|\,\nabla\boldsymbol{w}\,\\|_{\iota,h}.$

A combination of the above three inequalities, the discrete inf-sup condition (4.12) and the triangle inequalities gives (4.13).

Next, let $\boldsymbol{u}_{I}=I_{h}\boldsymbol{u}$ and let $p_{I}$ be the Clément interpolation [Clement:1975] of $p$ , respectively. The interpolation error (4.7) and the error estimates for the Clément interpolation give

|\!|\!|(\boldsymbol{u}-\boldsymbol{u}_{I},p-p_{I})|\!|\!|_{\iota,h}\leq Ch^{1/2}\Bigl{(}\|\,\boldsymbol{u}\,\|_{H^{3/2}}+\|\,p\,\|_{H^{1/2}}+\iota(\|\,\boldsymbol{u}\,\|_{H^{5/2}}+\|\,p\,\|_{H^{3/2}})\Bigr{)}\leq Ch^{1/2}\|\,\boldsymbol{f}\,\|_{L^{2}},

where we have used (2.29) and (2.30) in the last step.

Using the trace inequality (4.1), we bound the consistency error functional as

	$\displaystyle\iota^{2}\left\lvert\,\sum_{e\in\mathcal{E}_{h}}\int_{e}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot[\![\partial_{\boldsymbol{n}}\boldsymbol{w}]\!]\,\mathrm{d}\sigma(\boldsymbol{x})\,\right\rvert$	$\displaystyle\leq C\iota^{2}h^{1/2}(\\|\,\boldsymbol{u}\,\\|_{H^{2}}+\\|\,p\,\\|_{H^{1}})^{1/2}(\\|\,\boldsymbol{u}\,\\|_{H^{3}}+\\|\,p\,\\|_{H^{2}})^{1/2}\\|\,\nabla_{h}^{2}\boldsymbol{w}\,\\|_{L^{2}}$
		$\displaystyle\leq Ch^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}\\|\,\nabla\boldsymbol{w}\,\\|_{\iota,h},$

where we have used (2.28) in the last step.

Combining these inequalities, the discrete inf-sup condition (4.12) and the triangle inequalities gives (4.14). ∎

Corollary 4.6.

There exists $C$ independent of $\iota,\lambda$ and $h$ such that

(4.15)

|\!|\!|(\boldsymbol{u}_{0}-\boldsymbol{u}_{h},p_{0}-p_{h})|\!|\!|_{\iota,h}\leq C(\iota^{1/2}+h^{1/2})\|\,\boldsymbol{f}\,\|_{L^{2}},

where $\boldsymbol{u}_{0}$ is the solution of (2.8), and $p_{0}=\lambda\operatorname{div}\boldsymbol{u}_{0}$ .

Proof.

A combination of Theorem 2.9, Theorem 4.5, and the triangle inequality gives (4.15). ∎

5. Numerical Examples

In this part, we report the numerical performance for the proposed element of the lowest-order, i.e., $r=2$ . We test the accuracy and robustness of the element pair for the nearly incompressible materials. All examples are carried out on the nonuniform mesh. We are interested in the case when the Poisson’s ratio $\nu$ is close to $0.5$ and we report the relative errors $\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{h})\,\|_{\iota,h}/\|\nabla\boldsymbol{u}\|_{\iota}$ and the rates of convergence.

We let $\Omega=(0,1)^{2}$ , and set Young’s modulus $E=1$ . The Lamé constants are determined by

\lambda=\dfrac{E\nu}{(1+\nu)(1-2\nu)},\quad\mu=\dfrac{E}{2(1+\nu)}.

We set $\nu=0.3$ for the ordinary cases, hence $\lambda=0.5769$ , and $\mu=0.3846$ , and we set $\nu=0.4999$ for the nearly incompressible materials, hence $\lambda=\text{1.6664e3}$ , and $\mu=0.3334$ .

5.1. The first example

We test the performance of the element pair by solving a completely incompressible problem, which means $\operatorname{div}\boldsymbol{u}=0$ . Let $\boldsymbol{u}=(u_{1},u_{2})$ with

u_{1}=-\sin^{3}(\pi x)\sin(2\pi y)\sin(\pi y),\quad u_{2}=\sin(2\pi x)\sin(\pi x)\sin^{3}(\pi y).

Therefore $\operatorname{div}\boldsymbol{u}=0$ , and $\boldsymbol{f}$ is independent of $\lambda$ .

Table 2. Relative errors and convergence rates for the 1st example

$\iota\backslash h$	1/8	1/16	1/32	1/64
$\nu=0.3,\lambda=0.5769,\mu=0.3846$
1e+00	2.592e-01	1.333e-01	6.519e-02	3.246e-02
rate		0.96	1.03	1.01
1e-06	4.252e-02	1.159e-02	2.784e-03	6.918e-04
rate		1.88	2.06	2.01
$\nu=0.4999,\lambda=\text{1.6664e3},\mu=0.3334$
1e+00	2.592e-01	1.333e-01	6.519e-02	3.246e-02
rate		0.96	1.03	1.01
1e-06	4.252e-02	1.159e-02	2.784e-03	6.918e-04
rate		1.88	2.06	2.01

In view of Table 2, the optimal rates of convergence are observed with the completely incompressible media, which is consistent with the error bound (4.13).

5.2. The second example

This example is motivated by [Wihler:2006], which admits a singular solution. The exact solution $\boldsymbol{u}=(u_{1},u_{2})$ expressed in the polar coordinates as

u_{1}=u_{\rho}(\rho,\theta)\cos\theta-u_{\theta}(\rho,\theta)\sin\theta,\quad u_{2}=u_{\rho}(\rho,\theta)\sin\theta+u_{\theta}(\rho,\theta)\cos\theta,

where

	$\displaystyle u_{\rho}$	$\displaystyle=\dfrac{1}{2\mu}\rho^{\alpha}\Bigl{(}-(\alpha+1)\cos((\alpha+1)\theta)+(C_{2}-(\alpha+1))C_{1}\cos((\alpha-1)\theta)\Bigr{)},$
	$\displaystyle u_{\theta}$	$\displaystyle=\dfrac{1}{2\mu}\rho^{\alpha}\Bigl{(}(\alpha+1)\sin((\alpha+1)\theta)+(C_{2}+\alpha-1)C_{1}\sin((\alpha-1)\theta)\Bigr{)},$

and $\alpha=1.5,\omega=3\pi/4$ ,

C_{1}=-\dfrac{\cos((\alpha+1)\omega)}{\cos((\alpha-1)\omega)}\quad\text{and}\quad C_{2}=\dfrac{2(\lambda+2\mu)}{\lambda+\mu}.

It may be verified that $\boldsymbol{u}\in[H^{5/2-\varepsilon}(\Omega)]^{2}$ for a small number $\varepsilon>0$ . A direct calculation gives that $\boldsymbol{f}\equiv\boldsymbol{0}$ , while it is nearly incompressible because

\operatorname{div}\boldsymbol{u}=-\dfrac{3(1+\sqrt{2})}{\lambda+\mu}\rho^{1/2}\cos(\theta/2).

Table 3. Relative errors and convergence rates for the 2nd example

$\iota\backslash h$	1/8	1/16	1/32	1/64
$\nu=0.3,\lambda=0.5769,\mu=0.3846$
1e+00	1.062e-01	7.554e-02	5.355e-02	3.792e-02
rate		0.49	0.50	0.50
1e-06	2.809e-03	1.001e-03	3.549e-04	1.257e-04
rate		1.49	1.50	1.50
$\nu=0.4999,\lambda=\text{1.6664e3},\mu=0.3334$
1e+00	1.149e-01	8.200e-02	5.824e-02	4.135e-02
rate		0.49	0.49	0.49
1e-06	4.399e-03	1.567e-03	5.558e-04	1.968e-04
rate		1.49	1.50	1.50

It follows from Table 3 that the rates of convergence are sub-optimal. It is reasonable because the solution $\boldsymbol{u}$ is singular, which is similar to the results in [Wihler:2006]. The element pair is robust for the nearly incompressible materials.

5.3. The third example

In the last example, we test a problem with strong boundary layer effects. Such effects have been frequently observed in the stain elasticity model [Engel:2002, LiMingShi:2017, LiaoM:2019, LiMingWang:2021]. It is shown that the numerical solution converges to the solution of (2.8) when $\iota\ll h$ .

When $\iota\to 0$ , the boundary value problem (1.1) reduces to (2.8). Let $\boldsymbol{u}_{0}=(u_{1}^{0},u_{2}^{0})$ with

u_{1}^{0}=-\sin^{2}(\pi x)\sin(2\pi y),\quad u_{2}^{0}=\sin(2\pi x)\sin^{2}(\pi y)

be the solution of problem (2.8). The source term $\boldsymbol{f}$ is computed from (2.8). A direct calculation gives that $\operatorname{div}\boldsymbol{u}_{0}=0$ , and $\boldsymbol{f}$ is independent of $\lambda$ . The exact solution $\boldsymbol{u}$ for (1.1) is unknown, while it has strong boundary layer effects. In this case, we take $\iota\ll h$ , and report the relative error $\|\,\nabla(\boldsymbol{u}_{0}-\boldsymbol{u}_{h})\,\|_{\iota,h}/\|\nabla\boldsymbol{u}_{0}\|_{\iota}$ .

Table 4. Relative errors and convergence rates for the 3rd example

$\iota\backslash h$	1/8	1/16	1/32	1/64
$\nu=0.3,\lambda=0.5769,\mu=0.3846$
1e-04	1.311e-01	8.966e-02	6.299e-02	4.476e-02
rate		0.55	0.51	0.49
1e-06	1.311e-01	8.960e-02	6.283e-02	4.432e-02
rate		0.55	0.51	0.50
$\nu=0.4999,\lambda=\text{1.6664e3},\mu=0.3334$
1e-04	1.312e-01	8.968e-02	6.300e-02	4.476e-02
rate		0.55	0.51	0.49
1e-06	1.312e-01	8.963e-02	6.284e-02	4.432e-02
rate		0.55	0.51	0.50

It follows from Table 4 that the rate of convergence for the element pair changes to $1/2$ because of the boundary layer effects, which is consistent with the theoretical result. The element is still robust when the solution has strong boundary layer effects in the nearly incompressible limit.

	$\displaystyle\lambda\\|\,\nabla\operatorname{div}\boldsymbol{w}\,\\|_{L^{2}}^{2}$	$\displaystyle=\lambda(\nabla\operatorname{div}\boldsymbol{w},\nabla\operatorname{div}\boldsymbol{v}_{0})=A(\boldsymbol{w},\boldsymbol{v}_{0})-2\mu(\nabla\boldsymbol{\epsilon}(\boldsymbol{w}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))$
		$\displaystyle=(\boldsymbol{F},\boldsymbol{v}_{0})-2\mu(\nabla\boldsymbol{\epsilon}(\boldsymbol{w}),\nabla\boldsymbol{\epsilon}(\boldsymbol{v}_{0}))$
		$\displaystyle\leq\\|\,\boldsymbol{F}\,\\|_{H^{-2}}\\|\,\boldsymbol{v}_{0}\,\\|_{H^{2}}+2\mu\\|\,\nabla^{2}\boldsymbol{w}\,\\|_{L^{2}}\\|\,\nabla^{2}\boldsymbol{v}_{0}\,\\|_{L^{2}}$
		$\displaystyle\leq C\left(\\|\,\boldsymbol{F}\,\\|_{H^{-2}}+2\mu\\|\,\nabla^{2}\boldsymbol{w}\,\\|_{L^{2}}\right)\\|\,\nabla^{2}\boldsymbol{v}_{0}\,\\|_{L^{2}}$
		$\displaystyle\leq C\\|\,\boldsymbol{F}\,\\|_{H^{-2}}\\|\,\nabla\operatorname{div}\boldsymbol{w}\,\\|_{L^{2}}.$

	$\displaystyle 2\lambda\mu\iota^{2}\biggl{\|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}\boldsymbol{\epsilon}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{v}_{0})\mathrm{d}\sigma(\boldsymbol{x})\biggr{\|}\leq$	$\displaystyle C\Bigl{(}(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}})\\|\,p-p_{0}\,\\|_{\iota}$
		$\displaystyle\qquad+(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}})^{1/2}\\|\,p-p_{0}\,\\|_{\iota}^{3/2}\Bigr{)}$
	$\displaystyle\leq$	$\displaystyle\dfrac{1}{8}\\|\,p-p_{0}\,\\|_{\iota}^{2}+C\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota^{2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}\right).$

	$\displaystyle\iota^{2}\biggl{\|}\int_{\partial\Omega}\partial_{\boldsymbol{n}}pp_{0}\mathrm{d}\sigma(\boldsymbol{x})\biggr{\|}\leq$	$\displaystyle C\Bigl{(}\iota^{1/2}\\|\,p-p_{0}\,\\|_{\iota}\\|\,\boldsymbol{f}\,\\|_{L^{2}}+(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}+\iota^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}})^{3/2}\\|\,p-p_{0}\,\\|_{\iota}^{1/2}$
		$\displaystyle\qquad+\iota\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{3/2}\Bigr{)}$
		$\displaystyle\leq\dfrac{1}{8}\\|\,p-p_{0}\,\\|_{\iota}^{2}+C\left(\\|\,\nabla(\boldsymbol{u}-\boldsymbol{u}_{0})\,\\|_{\iota}^{2}+\iota\\|\,\boldsymbol{f}\,\\|_{L^{2}}^{2}\right).$

	$\displaystyle\iota^{2}\left\lvert\,\int_{\partial\Omega}(\partial_{\boldsymbol{n}}\boldsymbol{\sigma}\boldsymbol{n})\cdot\partial_{\boldsymbol{n}}\boldsymbol{u}_{0}\mathrm{d}\sigma(\boldsymbol{x})\,\right\rvert$	$\displaystyle\leq C\iota^{2}\left(\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{H^{1}}^{1/2}\\|\,\boldsymbol{\epsilon}(\boldsymbol{u}_{0})\,\\|_{H^{1}}+\\|\,\nabla\operatorname{div}\boldsymbol{u}\,\\|_{L^{2}}^{1/2}\\|\,\nabla\operatorname{div}\boldsymbol{u}\,\\|_{H^{1}}^{1/2}\\|\,p_{0}\,\\|_{H^{1}}\right)$
		$\displaystyle\leq C\iota^{2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\boldsymbol{u}\,\\|_{H^{3}}^{1/2}\left(\\|\,\boldsymbol{u}_{0}\,\\|_{H^{2}}+\\|\,p_{0}\,\\|_{H^{1}}\right)$
		$\displaystyle\leq C\iota^{2}\\|\,\nabla\boldsymbol{\epsilon}(\boldsymbol{u})\,\\|_{L^{2}}^{1/2}\\|\,\boldsymbol{u}\,\\|_{H^{3}}^{1/2}\\|\,\boldsymbol{f}\,\\|_{L^{2}}.$

	$\displaystyle\\|\,\nabla^{j}_{h}(\boldsymbol{v}-\varPi_{1}\boldsymbol{v})\,\\|_{L^{2}}\leq$	$\displaystyle\\|\,\nabla^{j}_{h}(\boldsymbol{v}-\varPi_{2}\boldsymbol{v})\,\\|_{L^{2}}+\\|\,\nabla^{j}_{h}\varPi_{0}(\boldsymbol{v}-\varPi_{2}\boldsymbol{v})\,\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{s-j}\\|\,\nabla^{s}\boldsymbol{v}\,\\|_{L^{2}}+Ch^{-j}\\|\,\boldsymbol{v}-\varPi_{2}\boldsymbol{v}\,\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{s-j}\\|\,\nabla^{s}\boldsymbol{v}\,\\|_{L^{2}}.$

Taylor-Hood like finite elements for nearly incompressible strain gradient elasticity problems

Abstract.

Key words and phrases:

2020 Mathematics Subject Classification:

1. Introduction

2. The Mixed Variational Formulation and Regularity Estimates

Example 2.1.

2.1. The mixed variational formulation

Lemma 2.2.

Proof.

Lemma 2.3.

Proof.

2.2. Regularity estimates

Lemma 2.4.

Proof.

Lemma 2.5.

Proof.

Lemma 2.6.

Proof.

Lemma 2.7.

Proof.

Lemma 2.8.

Proof.

Theorem 2.9.

Proof.

Corollary 2.10.

Proof.

3. A Family of Nonconforming Finite Elements

3.1. A family of finite elements

Lemma 3.1.

Proof.

Lemma 3.2.

Proof.

Remark 3.3.

3.2. Explicit representation for the bubble space

Lemma 3.4.

Proof.

Lemma 3.5.

Proof.

Example 3.6.

Lemma 3.7.

Proof.

Lemma 3.8.

Proof.

Example 3.9.

Remark 3.10.

4. The Mixed Finite Elements Approximation

4.1. The mixed finite element approximation

Lemma 4.1.

Proof.

Lemma 4.2.

Proof.

Theorem 4.3.

Proof.

4.2. Error estimates

Lemma 4.4.

Proof.

Theorem 4.5.

Proof.

Corollary 4.6.

Proof.

5. Numerical Examples

5.1. The first example

5.2. The second example

5.3. The third example

References