ANALYSIS OF A MODIFIED FIRST-ORDER ... - Semantic Scholar
ANALYSIS OF A MODIFIED FIRST-ORDER SYSTEM LEAST SQUARES METHOD FOR LINEAR ELASTICITY WITH IMPROVED MOMENTUM BALANCE GERHARD STARKE
∗,
† ¨ ¨ ALEXANDER SCHWARZ , AND JORG SCHRODER
Abstract. A modified first-order system least squares formulation for linear elasticity, obtained by adding the antisymmetric displacement gradient in the test space, is analyzed. This approach leads to surprisingly small momentum balance error compared to standard least squares approaches. It is shown that the modified least squares formulation is well-posed and its performance is illustrated by adaptive finite element computation based on using a closely related least squares functional as a posteriori error estimator. The results of our numerical computations show that, for the modified least squares approach, the momentum balance error converges at a much faster rate than the overall error. We prove that this is due to a strong connection of the stress approximation to that obtained from a mixed formulation based on the Hellinger-Reissner principle (with exact local momentum balance). The practical significance is that our proposed approach is almost momentum-conservative at a smaller number of degrees of freedom than mixed approximations with exact local momentum balance. Key words. First-order system least squares, linear elasticity, incompressible, momentum balance, a posteriori error estimator, Raviart-Thomas elements. AMS subject classifications. 65M60, 65M15
1. Introduction. In this paper, a modified first-order system least squares formulation for linear elasticity with improved momentum balance is analyzed. Its construction is obtained from the least squares approach in [11] by adding the antisymmetric displacement gradient in the test space. This formulation was introduced in [15] from an engineering perspective and shown to give much better results in bending dominated situations. The resulting variational formulation is based on a nonsymmetric bilinear form. The ellipticity properties of this bilinear form are shown to be essentially the same as for the symmetric bilinear form associated with the “true” least squares formulation from [11]. The numerical results in [15] show a significantly improved momentum balance for our nonsymmetric formulation compared to the standard least squares method. This seems to give an explanation for the improved results, especially for bending dominated problems. In this paper, we are able to show that the reduction rate for the momentum balance error is indeed of higher order than the overall approximation error. The improvement of the reduction rate depends on the regularity of the problem. A major step in our analysis is the derivation of a close connection between the stress approximation of our method and that obtained by a mixed formulation using the Arnold-Winther element ([3], generalized to three dimensions in [1]). Such a connection was established in [7] for the least squares formulation of second-order elliptic boundary value problem and its related saddle point and primal formulations. Least squares finite element methods became very popular in recent years in many different application areas, see [5]. An important property of mixed formulations based on least squares functionals is the fact that its local evaluation provides a local a posteriori error estimate (see e.g. [4] or, in the context of linear elasticity ∗ Institut f¨ ur Angewandte Mathematik, Leibniz Universit¨ at Hannover, Welfengarten 1, 30167 Hannover, Germany ([email protected]). † Institut f¨ ur Mechanik, Fakult¨ at f¨ ur Ingenieurwissenschaften, Universit¨ at Duisburg-Essen, Universit¨ atsstr. 15, 45117 Essen, Germany ({alexander.schwarz,j.schroeder}@uni-due.de).
1
2
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
[9]). Most importantly, such approaches simultaneously construct approximations to all variables of interest (like displacements and stresses) in appropriate spaces. Alternatively, mixed methods based on saddle point formulations may be used for this purpose. For this approach, however, suitable combinations of finite element spaces are restricted by the inf-sup conditions. Families of finite element spaces for such mixed formulations of saddle point type may be found e.g. in [12]. Those families of finite element spaces are based on Raviart-Thomas elements (as a generalization of PEERS, see [2]) or Brezzi-Douglas-Marini elements (as a generalization of BDMS, see [16]). A posteriori error estimators were constructed in [13] for these saddle point approaches based on the Helmholtz decomposition. Of course, the Arnold-Winther elements from [3] which are used in this contribution as a theoretical tool could in principle also be chosen. However, these finite element spaces involve an even higher number of unknowns, in general. The outline of this paper is as follows. In the following section our modified least squares formulation for the linear elasticity model is introduced. Section 3 contains the ellipticity result for our new variational approach. We proceed with the finite element formulation and the results of our numerical computations on uniformly and adaptively refined triangulations in section 4. Finally, the strong connection between our new formulation and the mixed formulation of saddle point type based on the Arnold-Winther elements as well as the standard displacement formulation is established in section 5. 2. The Modified First-Order System Least Squares Formulation. We start from the equations of linear elasticity in the form div σ = 0 , σ − Cε(u) = 0 ,
(2.1)
where ε(u) = (∇u + ∇uT )/2 denotes the linear strain tensor and C describes the linear material law given by ν E (tr ε(u)) I ε(u) + Cε(u) = 1+ν 1 − 2ν with the elasticity modulus E and Poisson ratio ν. The boundary of Ω ⊂ IRd is divided into two parts as ∂Ω = ΓD ∪ ΓN . We consider boundary conditions of the form u = 0 on ΓD , σ · n = g on ΓN and define 1 HD (Ω) = {q ∈ H 1 (Ω) : q = 0 on ΓD }
HN (div, Ω) = {w ∈ H(div, Ω) : w · n = 0 on ΓN } . Throughout this paper, the norm on L2 (Ω) (or L2 (Ω)d , L2 (Ω)d×d , respectively) will be abbreviated by k · k and the corresponding inner product by (·, ·). We assume that ΓD is a subset of ∂Ω with positive measure (length if d = 2, area in the three-dimensional case) such that Korn’s inequality is valid in the form 1 kvk2 + k∇vk2 ≤ CK kC 1/2 ε(v)k2 for all v ∈ HD (Ω)d
(2.2)
with a constant CK > 0 (cf. [6, Section VI.3]). Note that CK is independent of ν due to the fact that kC 1/2 ε(v)k2 =
E ν E kε(u)k2 + ktr ε(u))k2 ≥ kε(u)k2 . 1+ν 1 − 2ν 1+ν
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
3
We will frequently use the decomposition of an arbitrary matrix-valued function τ ∈ L2 (Ω)d×d into its symmetric and antisymmetric part, τ = sy τ + as τ with sy τ =
τ + τT τ − τT , as τ = . 2 2
Obviously, (sy τ , as τ ) = 0 which implies kτ k2 = ksy τ k2 + kas τ k2 . If σ N ∈ H(div, Ω)d is chosen such that σ N · n = g on ΓN , then the solution of (2.1) may be obtained from minimizing the least squares functional G(σ, u) = kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2
(2.3)
1 among all σ ∈ σ N + HN (div, Ω)d and u ∈ HD (Ω)d . This is exactly the least squares functional studied in [9]. One may also extend this functional by the redundant equation
as (C −1/2 σ) = 0 (cf. [10]) which leads to F(σ, u) = kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2 + kas (C −1/2 σ)k2 .
(2.4)
The minimization of (2.4) is equivalent to the variational problem of finding σ ∈ 1 σ N + HN (div, Ω)d and u ∈ HD (Ω)d such that (div σ, div τ ) + (C −1 σ − ε(u), τ ) + (as (C −1 σ), as τ ) = 0 , (σ − C ε(u), ε(v)) = 0
(2.5)
1 (Ω)d . The additional term associated with holds for all τ ∈ HN (div, Ω)d and v ∈ HD the antisymmetric stress part is not necessary in the above least squares formulation but will be required in the following modification. Replacing ε(v) by ∇v in the displacement variation of (2.5) leads to our modified least squares method. It consists 1 in finding σ ∈ σ N + HN (div, Ω)d and u ∈ HD (Ω)d such that
(div σ, div τ ) + (C −1 σ − ε(u), τ ) + (as (C −1 σ), as τ ) = 0 , (σ − Cε(u), ∇v) = 0
(2.6)
1 holds for all τ ∈ HN (div, Ω)d and v ∈ HD (Ω)d .
3. Ellipticity. We consider the nonsymmetric bilinear form B(σ, u; τ , v) = (div σ, div τ ) + (C −1 σ − ε(u), τ − C ∇v) + (as (C −1 σ), as τ ) (3.1) associated with the modified least squares formulation (2.6). We will show that this bilinear form is coercive and continuous with respect to the product space HN (div, Ω)d × 1 HD (Ω)d equipped with the norms 1/2 kτ kdiv,C −1 = kdiv τ k2 + kC −1/2 τ k2 , 1/2 kvk1,C = kC 1/2 ∇vk2 + kvk2 .
(3.2)
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
4
For these estimates we rely on the assumption that E is on the order of one. Note that this can always be achieved by a suitable rescaling of the variables. Coercivity and continuity of the bilinear form is then achieved uniformly in the incompressible limit, i.e. as ν → 1/2. This is the statement of the following theorem. Theorem 3.1. There exist positive constants α and β which are independent of ν such that B(τ , v; τ , v) ≥ α kτ k2div,C −1 + kvk21,C (3.3) 1 holds for all (τ , v) ∈ HN (div, Ω)d × HD (Ω)d and
1/2 1/2 kτ k2div,C −1 + kvk21,C B(σ, u; τ , v) ≤ β kσk2div,C −1 + kuk21,C
(3.4)
1 holds for all (τ , v), (σ, u) ∈ HN (div, Ω)d × HD (Ω)d . Proof. The proof of (3.4) is rather straightforward. The Cauchy-Schwarz inequality immediately implies 1/2 B(σ, u; τ , v) ≤ kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2 + kas (C −1/2 σ)k2 (3.5) 1/2 kdiv τ k2 + kC −1/2 τ − C 1/2 ∇vk2 + kas (C −1/2 τ )k2 .
Using the orthogonal splitting of an arbitrary matrix into its symmetric and antisymmetric part, we obtain kC 1/2 ∇vk2 = (C ∇v, ∇v) = (sy (C ∇v), ε(v)) + (as (C ∇v), as ∇v) = (C ε(v), ε(v)) + (C as ∇v, as ∇v) = kC 1/2 ε(v)k2 + kC 1/2 as ∇vk2 , which holds due to E ν sy ∇v + (tr ∇v) I 1+ν 1 − 2ν E ν = ε(v) + (tr ε(v)) I = Cε(v) , 1+ν 1 − 2ν E as (C ∇v) = as ∇v = C as ∇v . 1+ν
sy (C ∇v) =
In an analogous way, one obtains kC −1/2 τ k2 = kas (C −1/2 τ )k2 + ksy (C −1/2 τ )k2 . These identities imply kC 1/2 ε(v)k ≤ kC 1/2 ∇vk , kas C −1/2 τ k ≤ kC −1/2 τ k ,
(3.6)
respectively. Using (3.6), (3.5) can be further estimated as 1/2 B(σ, u; τ , v) ≤ kdiv σk2 + 2kC −1/2 σk2 + 2kC 1/2 ∇uk2 + kC −1/2 σk2 1/2 kdivτ k2 + 2kC −1/2 τ k2 + 2kC 1/2 ∇vk2 + kC −1/2 τ k2 .
(3.7)
5
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
This means that (3.4) holds with β = 3. Naturally, the coercivity estimate (3.3) is harder to establish. We start from B(τ , v; τ , v) = kdiv τ k2 + (C −1 τ − ε(v), τ − C ∇v) + kas (C −1/2 τ )k2 and note that all three terms on the right hand side are nonnegative. This is trivial for the first and third term. The second term constitutes a quadratic functional which 1 clearly attains its minimum for (τ ∗ , v∗ ) ∈ HN (div, Ω)d × HD (Ω)d satisfying 2(C −1 τ ∗ , τ ) − (ε(v∗ ) + ∇v∗ , τ ) = 0 for all τ ∈ HN (div, Ω)d , 1 2(C ε(v∗ ), ε(v)) − (τ ∗ , ε(v) + ∇v) = 0 for all v ∈ HD (Ω)d .
In particular, 2(C −1 τ ∗ , τ ∗ ) − (ε(v∗ ) + ∇v∗ , τ ∗ ) = 0 , 2(C ε(v∗ ), ε(v∗ )) − (τ ∗ , ε(v∗ ) + ∇v∗ ) = 0 which implies (C −1 τ ∗ − ε(v∗ ), τ ∗ − C ∇v∗ ) = (C −1 τ ∗ , τ ∗ ) − (τ ∗ , ∇v∗ + ε(v∗ )) + (C ε(v∗ ), v∗ ) = 0 . With a constant C1 ≥ 1 which is still free to be chosen appropriately, we therefore have 1 kdiv τ k2 + (C −1 τ − ε(v), τ − C ∇v) + C1 kas (C −1/2 τ )k2 B(τ , v; τ , v) ≥ C1 1 = kdiv τ k2 + (C −1 τ − ε(v), τ − C ε(v)) C1 (3.8) −(C −1 τ − ε(v), as (C ∇v)) + C1 kas C −1/2 τ k2 1 = kdiv τ k2 + kC −1/2 τ − C 1/2 ε(v)k2 C1 −(as (C −1/2 τ ), as (C 1/2 ∇v)) + C1 kas (C −1/2 τ )k2 . We will rely on the coercivity proof for the symmetric bilinear form from [9, Theorem 2.1] which states that, with a positive constant C0 (independently of ν), G(τ , v) ≥ C0 kdiv τ k2 + kC −1/2 τ k2 + kC 1/2 ε(v)k2 (3.9) 1 holds for all (τ , v) ∈ HN (div, Ω)d × HD (Ω)d . From Korn’s inequality (2.2) we deduce
E E kas ∇vk2 ≤ k∇vk2 1+ν 1+ν CK E CK E ν ≤ kε(v)k2 ≤ kε(v)k2 + kdiv vk2 1+ν 1+ν 1 − 2ν
kas (C 1/2 ∇v)k2 = (as (C ∇v), as ∇v) =
(3.10)
= CK (C ε(v), ε(v)) = CK kC 1/2 ε(v)k2 . Combining (3.10) with (3.9) leads to kas (C 1/2 ∇v)k2 ≤
CK G(τ , v) . C0
(3.11)
6
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
Therefore, (3.8) can be further bounded from below as follows: 1 B(τ , v; τ , v) ≥ kdivτ k2 + kC −1/2 τ − C 1/2 ε(v)k2 + C1 kas C −1/2 τ k2 C1 −(as (C −1/2 τ ), as (C 1/2 ∇v)) 1 kdiv τ k2 + kC −1/2 τ − C 1/2 ε(v)k2 + C1 kas (C −1/2 τ )k2 ≥ C1 1 1/2 2 −1/2 2 kas (C ∇v)k −C1 kas (C τ )k − 4C1 1 1 1/2 2 = G(τ , v) − kas (C ∇v)k C1 4C1 1 CK ≥ 1− G(τ , v) . C1 4C0 C1
(3.12)
If we choose C1 = max{1, CK /(2C0 )}, we have CK /(4C0 C1 ) ≤ 1/2 and therefore B(τ , v; τ , v) ≥
1 G(τ , v) . 2C1
Combined with (3.9) and (2.2), this leads to C0 1 2 −1/2 2 2 1/2 2 B(τ , v; τ , v) ≥ kdiv τ k + kC τk + kvk + kC ∇vk 2C1 CK
(3.13)
(3.14)
which finishes the proof. The constants α and β in Theorem 1 are independent of ν assuming that E is on the order of 1. This means that our coercivity and continuity estimates hold uniformly in the incompressible limit as ν → 1/2. Coercivity and continuity could also be established with respect to unscaled norms 1/2 k · kdiv := k · kdiv,I = kdiv ( · ) k2 + k · k2 , 1/2 k · k1 := k · k1,I = k∇ ( · ) k2 + k · k2 . However, these estimates would not be uniform in the incompressible limit. The fact that (3.3) and (3.4) hold ensures that our modified least squares formulation is well-posed, i.e. the variational problem (2.6) has a unique solution which depends in a continuous way on the boundary data. 4. Finite Element Approximation. For the finite element approximation of the variational problem (2.6), we choose finite-dimensional subspaces Σh ⊂ Σ := 1 HN (div, Ω)d and Vh ⊂ V := HD (Ω)d . The discrete variational problem consists in finding σ h = σ N + Σh and uh ∈ Vh such that (div σ h , div τ h ) + (C −1 σ h − ε(uh ), τ h − C ∇vh ) + (as (C −1 σ h ), as τ h ) = 0
(4.1)
holds for all τ h ∈ Σh and vh ∈ Vh . Combined with (2.6), this implies that (div (σ − σ h ), div τ h ) + (C −1 (σ − σ h ) − ε(u − uh ), τ h − C ∇vh ) + (as (C −1 (σ − σ h )), as τ h ) = 0 is satisfied for all τ h ∈ Σh and vh ∈ Vh .
(4.2)
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
7
In our computational results, next-to-lowest order Raviart-Thomas elements (RT1 ) for the stress space Σh are combined with conforming quadratic elements (P2 ) for the displacement space Vh . Example 1. We consider Cook’s membrane problem with corners (0, 0), (48, 44), (48, 60) and (0, 44) which is clamped at the left boundary segment (see figure 4.1). The rest of the boundary is subjected to surface traction forces which are zero on the upper and lower boundary segments and pointing upwards on the right boundary segment. We consider an almost incompressible material with Poisson ratio ν = 0.49.
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 8 32 128 512 2048 8192
dim Σh 40 144 544 2112 8320 33024
dim Vh 72 304 1248 5056 20352 81664
F(σ h , uh ) 12.729 5.402 2.642 1.206 0.617 0.421
kdiv σ h k2 2.640 2.125 1.240 0.566 0.387 0.334
Table 4.1 Symmetric least squares formulation (2.5): RT1 × P2 , uniform refinement
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 8 32 128 512 2048 8192
dim Σh 40 144 544 2112 8320 33024
dim Vh 72 304 1248 5056 20352 81664
F(σ h , uh ) 19.145 9.421 4.468 1.631 0.753 0.527
kdiv σ h k2 4.527 2.336 0.901 0.499 0.374 0.289
Table 4.2 Nonsymmetric least squares formulation (2.6): RT1 × P2 , uniform refinement
A comparison of the results in tables 4.1 and 4.2 shows that, on the same triangulations resulting from uniform refinement, the nonsymmetric formulation produces a faster decrease of the momentum balance error. Note that the overall least squares functional must be smaller for the symmetric formulation (on identical triangulations) since it is actually minimized. However, the higher reduction rate of the momentum balance error is not very pronounced in this example. We also compare these two formulations on a sequence of adaptively refined triangulations based on the a posteriori error estimator provided by the local evaluation of the functional F(σ h , uh ). A fixed proportion of elements with largest local contribution to the functional is refined in our adaptive approach (about 15 % of elements in our computations). The results in tables 4.3 and 4.4 show that on adaptively refined triangulations, the difference between the two formulations is much more striking. The reduction rate of the momentum balance term in the functional is significantly faster in the nonsymmetric formulation. Surprisingly, even the least squares functional on the finest refinement level is smaller for the nonsymmetric formulation. This means that these refined triangulations are better suited to our problem. The adaptively refined triangulations as well as the deformed domain are shown in figure 4.1. The most striking difference is that the nonsymmetric formulation (on
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
8
l=0 l=1 l=2 l=3 l=4 l=5 l=6 l=7 l=8 l=9 l = 10
# elem. 8 19 43 90 184 372 732 1359 2440 4351 7568
dim Σh 40 90 194 396 796 1572 3044 5586 9960 17634 30600
dim Vh 72 176 408 864 1780 3636 7204 13440 24200 43280 75352
F(σ h , uh ) 12.729 5.476 2.584 1.627 1.061 0.748 0.519 0.413 0.341 0.267 0.185
kdiv σ h k2 2.640 2.048 1.257 0.975 0.718 0.543 0.394 0.323 0.253 0.162 0.082
Table 4.3 Symmetric least squares formulation (2.5): RT1 × P2 , adaptive refinement
l=0 l=1 l=2 l=3 l=4 l=5 l=6 l=7 l=8 l=9 l = 10
# elem. 8 19 41 87 190 370 710 1346 2464 4415 7774
dim Σh 40 90 186 382 812 1560 2948 5528 10044 17902 31408
dim Vh 72 176 388 836 1848 3620 6992 13316 24452 43908 77428
F(σ h , uh ) 19.145 10.773 4.595 2.558 1.277 0.796 0.592 0.476 0.360 0.230 0.129
kdiv σ h k2 4.527 3.361 1.874 1.010 0.548 0.379 0.269 0.151 0.065 0.020 0.005
Table 4.4 Nonsymmetric least squares formulation (2.6): RT1 × P2 , adaptive refinement
the right) puts much more emphasis at refining near the upper left corner. Example 2. We consider the more regular Cook’s membrane problem with corners (0, 0), (48, 14), (48, 30) and (0, 44) where the singularities are much less severe. This is a good illustration of our analysis in the next section which will rely on the regularity of the underlying boundary value problem. It will indicate that the decrease of the momentum balance error relative to the overall error will be more pronounced the more regular the problem is. We repeat our computations for this more regular modification of Cook’s membrane on a sequence of uniformly refined triangulations. The results in table 4.6 for the nonsymmetric formulation clearly show a faster reduction of the momentum balance error which the numbers in table 4.5 for the symmetric formulation do not. 5. A Strong Connection To Mixed and Displacement Formulations. In [7], Brandts, Chen and Yang have established a strong connection for second-order elliptic problems between the least squares formulation and corresponding mixed and standard Galerkin approaches. Such a strong connection also holds for the linear elasticity problem between our modified variational formulation (2.6) and an appropriate mixed and the standard displacement-based approach. This result is the basis for
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
9
Fig. 4.1. Triangulation of deformed Cook’s membrane after seven adaptive refinement steps: symmetric (left) vs. nonsymmetric (right) formulation
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 11 44 176 704 2816 11264
dim Σh 54 196 744 2896 11424 45376
dim Vh 100 420 1720 6960 28000 112320
F(σ h , uh ) 4.575 1.890 0.892 0.511 0.338 0.204
kdiv σ h k2 1.491 0.789 0.467 0.330 0.197 0.079
Table 4.5 Symmetric least squares formulation (2.5): RT1 × P2 , uniform refinement
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 11 44 176 704 2816 11264
dim Σh 54 196 744 2896 11424 45376
dim Vh 100 420 1720 6960 28000 112320
F(σ h , uh ) 7.662 2.551 1.151 0.724 0.487 0.264
kdiv σ h k2 3.149 0.916 0.397 0.180 0.056 0.011
Table 4.6 Nonsymmetric least squares formulation (2.6): RT1 × P2 , uniform refinement
our proof of the improved momentum balance observed in the computational results shown in Section 4. In our analysis, we will use the finite element spaces introduced by Arnold and
10
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
Winther in [3] which are defined as follows for k ≥ 1: 2×2 ΣAW = {τ ∈ Σ : as τ = 0 , τ |T ∈ Pk+2 , div τ |T ∈ Pk2 } , h
(5.1)
where Pk denotes the space of polynomials of degree k (in two dimensions). Important for our consideration is the fact that div ΣAW coincides with the space of piecewise h polynomials of degree k (without any inter-element continuity conditions). We restrict our analysis to two dimensional problems. For some fixed integer k ≥ 1, let ΣAW ⊆ Σ denote the Arnold-Winther space of h degree k introduced above. Moreover, let Vh ⊂ V be the standard H 1 (Ω)-conforming finite element space of piecewise polynomials of degree k + 1. Finally, denote by Xh ⊂ L2 (Ω)2 the piecewise polynomial spaces of degree k without any continuity requirements. Let us denote the finite element approximation based on the variational formuAW AW AW N lation (4.1) with Σh = ΣAW by (σ AW + h h , uh ). Note that (σ h , uh ) ∈ (σ AW Σh ) × Vh satisfies as σ h = 0 and therefore also minimizes the least squares functional G(σ h , uh ) with respect to ΣAW h ) × Vh . The first step of our analysis consists AW in establishing a close connection of (σ AW h , uh ) to the following mixed and standard displacement approaches. The mixed formulation that we consider is based on Hellinger-Reissner principle and uses the finite element space by Arnold and Winther [3]. The mixed finite element m AW N approximation (σ m h , ξ h ) ∈ (σ + Σh ) × Xh with these spaces satisfy m (C −1 (σ − σ m h ), τ h ) + (u − ξ h , div τ h ) = 0
(5.2)
for all τ h ∈ ΣAW and h (div(σ − σ m h ), η h ) = 0 for all η h ∈ Xh .
(5.3)
For the combination ΣAW × Xh introduced above (Arnold-Winther space combined h with piecewise polynomials of degree k) the first constraint in (5.3) implies that 2 div σ m h = 0 holds exactly in L (Ω). Finally, the standard displacement formulation udh ∈ Vh satisfies (Cε(u − udh ), ε(vh )) = 0
(5.4)
for all vh ∈ Vh . We are now ready to prove the generalization of the Brandts-Chen-Yang formula (cf. [7, Lemma 3.2]) to our modified least squares formulation (2.6). Lemma 5.1. For some integer k ≥ 1, let ΣAW denote the Arnold-Winther space h of degree k and let Vh denote the standard continuous elements of degree k + 1. For AW AW N the finite element approximation (σ AW h , uh ) ∈ (σ + Σh ) × Vh based on (4.1), the AW m N mixed approximation σ h ∈ σ +Σh based on (5.2) and (5.3), and the displacement approximation udh ∈ Vh based on (5.4), the following relation holds: AW AW B(σ AW − σm − udh ; σ AW − σm − udh ) h h , uh h h , uh AW d AW = (C −1 (σ − σ m − σm − σm h ), σ h h ) − (ε(u − uh ), σ h h ).
(5.5)
11
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
Proof. The definition of the bilinear form (3.1) leads to AW AW B(σ AW − σm − udh ; σ AW − σm − udh ) h h , uh h h , uh 2 = kdiv (σ AW − σm h h )k AW AW + (C −1 (σ AW − σm − udh ), σ AW − σm − udh )) h h ) − ε(uh h h − C ∇(uh AW + (as (σ AW − σm − σm h h ), as (σ h h ))
= (div (σ AW − σ), div (σ AW − σm h h h ) AW + (C −1 (σ AW − σ) − ε(uAW − u), σ AW − σm − udh )) h h h h − C ∇(uh
(5.6)
+ (as (σ AW − σ), as (σ AW − σm h h h )) d AW AW + (C −1 (σ − σ m − σm − udh )) h ) − ε(u − uh ), σ h h − C ∇(uh AW + (as (σ − σ m − σm h ), as (σ h h )) AW m AW = (C −1 (σ − σ m − σm − udh )) h ), σ h h ) − (σ − σ h , ∇(uh d AW − (ε(u − udh ), σ AW − σm − udh )) h h ) + (C ε(u − uh ), ε(uh AW + (as (σ − σ m − σm h ), as (σ h h )) ,
AW m in the last equality (note that σ AW where we used (4.2) with Σh = ΣAW h h −σ h ∈ Σh and uAW − udh ∈ Vh ). The second and fourth terms in (5.6) vanish due to (5.3) and h (5.4), respectively. The last term vanishes since symmetry is imposed in the stress space ΣAW h . This leads us to the desired generalization of the Brandts-Chen-Yang formula (5.5). In the sequel, we will use discrete inverses of the divergence operator with respect to the spaces ΣRT ⊆ Σ and ΣAW ⊆ Σ (cf. [14, Sect. 2]). To this end, for any h h AW to be the solution of the saddle point problem ∈ Σ ξ h ∈ Xh , we may define ΞAW h h
(C −1 ΞAW h , τ h ) − (ξ h , div τ h ) = 0
for all τ h ∈ ΣAW , h
(div ΞAW h , η h ) = (ξ h , η h )
for all η h ∈ Xh .
(5.7)
RT Similarly, we may define ΞRT h ∈ Σh to be the solution of −1 RT Ξh ), as τ h ) − (ξ h , div τ h ) = 0 for all τ h ∈ ΣRT (C −1 ΞRT h , h , τ h ) + (as(C
(div ΞRT h , η h ) = (ξ h , η h ) for all η h ∈ Xh .
(5.8)
If we define subspaces RT/AW,0
Σh
RT/AW
= {τ h ∈ Σh
: div τ h = 0}
(5.9)
and the inner product h · , · iC −1 ,as on Σ given by h σ , τ iC −1 ,as = (C −1 σ, τ ) + (as(C −1 σ), as τ ) ,
(5.10)
AW,0 we obtain hΞAW and, similarly, hΞRT h , τ h iC −1 ,as = 0 for all τ h ∈ Σh h , τ h iC −1 ,as = 0 RT,0 for all τ h ∈ Σh . Furthermore, we define orthogonal subspaces RT/AW,⊥
Σh
RT/AW
= {τ h ∈ Σh
RT/AW,0
: hτ h , ω h iC −1 ,as = 0 for all ω h ∈ Σh
}.
(5.11)
The construction in (5.7) and (5.8) provides discrete inverses of the divergence operRT/AW,⊥ ator mapping Σh , respectively, onto Xh . We therefore use the notation −1 AW ΞRT = div−1 h = divRT ξ h and Ξh AW ξ h .
(5.12)
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
12
For any ξ ∈ L2 (Ω)2 , we may proceed in any similar way as in (5.7) and (5.8) with respect to the entire Sobolev space Σ = HN (div, Ω)2 . To this end, let Ξ ∈ Σ be the solution of (C −1 Ξ, τ ) + (as(C −1 Ξ), as τ ) − (ξ, div τ ) = 0 (div Ξ, η) = (ξ, η)
for all τ ∈ Σ , for all η ∈ L2 (Ω)2 .
(5.13)
This defines the inverse of the divergence operator mapping the orthogonal subspace Σ⊥ of Σ with respect to Σ0 = {τ ∈ Σ : div τ = 0}
(5.14)
onto L2 (Ω)2 , Ξ = div−1 ξ. From now on, we will use the notation ιh . κh to denote that ιh ≤ Cκh holds with some constant C independently of h. Some fundamental properties of the operators introduced in (5.7) and (5.8) are summarized in the following lemma. The statements of Lemma 5.2 are well-known and, in the case of the Arnold-Winther spaces, also used in [14]. We state them explicitly here for the sake of better readability. Lemma 5.2. Under the assumptions that the problem (5.13) satisfies the regularity assumption Ξ ∈ H 1 (Ω)2×2 with kΞk1,Ω . kξk ,
(5.15)
the following holds for the discrete inverses of the divergence operator defined in (5.7) and (5.8): −1 kdiv−1 ξ h − div−1 RT ξ h k . h kξ h k , kdivRT ξ h k . kξ h k
(5.16)
for all ξ h ∈ Xh and −1 kdiv−1 ξ h − div−1 AW ξ h k . h kξ h k , kdivAW ξ h k . kξ h k
(5.17)
for all ξ h ∈ Xh . Proof. The regularity assumption (5.15) implies kΞk1,Ω . kdiv Ξk .
(5.18)
The standard approximation properties of Raviart-Thomas spaces (see [8, Prop. 3.9]) lead to kΞ − ΞRT h k . h kΞk1,Ω which, combined with (5.18), proves (5.16). Finally, the first approximation property stated in [3, Theorem 6.1] implies kΞ − ΞAW h k . h kΞk1,Ω which, combined with (5.18), leads to (5.17). The next step consists in showing that the momentum balance error of the formulation (4.1) using Arnold-Winther elements converges at a faster rate than the error norm measured by the least squares functional.
13
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
m AW d Theorem 5.3. Let σ AW h , σ h , uh , uh be as in Lemma 5.1 and assume that (5.15) is satisfied. Furthermore, assume that the boundary value problem associated with (2.1) is H 2 -regular in the sense that, for any g ∈ H 1/2 (ΓN )2 ,
u ∈ H 2 (Ω)2 with kuk2,Ω . kgk1/2,ΓN
(5.19)
holds. Then, 1/2 −1/2 2 1/2 kdiv σ AW (σ − σ m (ε(u) − ε(udh ))k2 h k . h kC h )k + kC
(5.20)
and AW AW 1/2 kdiv σ AW . h k . h F(σ h , uh )
(5.21)
Proof. From (5.5) and (3.3), we obtain (note that div σ m h = 0) 2 AW 2 kdiv σ AW − σm h k = kdiv (σ h h )k AW d AW . (C −1 (σ − σ m − σm − σm h ), σ h h ) − (ε(u − uh ), σ h h ).
(5.22)
For the first term in (5.22), (5.2) leads to m AW AW (C −1 (σ − σ m − σm − σm h ), σ h h ) = (u − ξ h , div (σ h h )) m AW AW = (u − ξ m h , div σ h ) = (u − ξ h , div Ξh )
= (C
−1
(σ −
AW σm h ), Ξh )
(5.23)
,
= is the solution of (5.7) with ξ h = div σ AW (note that div ΞAW ∈ ΣAW where ΞAW h h h h AW div σ h holds). We also have −1 −1 (C −1 (σ − σ m (σ − σ m (σ − σ m h ), Ξ) = (C h ), Ξ) + (as (C h )), as Ξ) = 0 ,
(5.24)
if Ξ ∈ Σ is the solution of (5.13) with ξ = div σ AW (note that div (σ − σ m h h ) = 0 and m as (σ − σ h ) = 0 holds). The combination of (5.23) and (5.24) leads to AW AW −1 − Ξ) (C −1 (σ − σ m − σm (σ − σ m h ), σ h h ) = (C h ), Ξh −1/2 ≤ kC −1/2 (σ − σ m (Ξ − ΞAW h )kkC h )k
. h kC
−1/2
(σ −
σm h )k
kdiv
σ AW h k
(5.25)
,
AW where kC −1/2 (Ξ − ΞAW h )k . hkdiv σ h k is used which holds due to Lemma 5.2. For the second term in (5.22), integration by parts gives d AW d AW −(ε(u − udh ), σ AW − σm − σm h h ) = (u − uh , div (σ h h ) = (u − uh , div σ h ) ,
where as (σ AW − σm h h ) = 0 was used. This leads to d AW −(ε(u − udh ), σ AW − σm h h ) ≤ ku − uh k kdiv σ h k
. hkC 1/2 ε(u − udh )k kdiv σ AW h k,
(5.26)
where the estimate ku−udh k . hkC 1/2 ε(u−udh )k was used which holds due to the classical Aubin-Nitsche duality argument and (2.2) under our assumption of H 2 regularity for u. The combination of (5.25) and (5.26) proves (5.20).
14
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
The estimate (5.21) is obtained from (5.20) in the following way: The relation (5.2) implies 2 −1/2 kC −1/2 (σ − σ m (σ − τ h )k2 h )k ≤ kC
for all τ h ∈ ΣAW,0 and therefore h 2 −1/2 AW 2 kC −1/2 (σ − σ m (σ − (σ AW − div−1 h )k ≤ kC h AW div σ h )k −1 2 AW 2 ≤ 2 kC −1/2 (σ − σ AW )k + kdiv div σ )k h h AW 2 AW 2 . 2 kC −1/2 (σ − σ AW h )k + kdiv σ h k 2 AW 2 , = 2 kC −1/2 (σ − σ AW )k + kdiv (σ − σ )k h h
(5.27)
where Lemma 5.2 is used again. Similarly, the relation (5.4) leads to kC 1/2 (ε(u) − ε(udh ))k2 ≤ kC 1/2 (ε(u) − ε(vh ))k2 for all vh ∈ Vh which implies 2 kC 1/2 (ε(u) − ε(udh ))k2 ≤ kC 1/2 (ε(u) − ε(uAW h ))k .
(5.28)
From (5.27), (5.28), together with (3.9), 2 1/2 AW AW AW kC −1/2 (σ − σ m (ε(u) − ε(udh ))k2 . G(σ AW h )k + kC h , uh ) = F(σ h , uh ) (5.29)
is obtained which finishes the proof of (5.21). RT AW AW N RT N Theorem 5.4. Let (σ RT h , uh ) ∈ (σ + Σh ) × Vh and (σ h , uh ) ∈ (σ + AW Σh ) × Vh be the approximations of our modified least squares method (4.1) in the respective finite element spaces. Then, under the assumption that (5.15) holds, AW RT RT 1/2 kdiv σ RT . h k . kdiv σ h k + h F(σ h , uh )
(5.30)
Proof. We start with the observation that divσ AW as well as divσ RT h h are contained in Xh . Let AW ⊥ div−1 = PAW σ AW , h AW div σ h RT ⊥ RT div−1 RT div σ h = PRT σ h
(5.31)
be the orthogonal projections with respect to h · , · iC −1 ,as onto ΣAW,⊥ and ΣRT,⊥ , h h AW,⊥ −1 respectively. Plugging the special test functions τ h = divAW ξ h ∈ Σh with ξh ∈ Xh into (4.1) leads to −1 −1 AW AW (div σ AW h , ξ h ) + hσ h , divAW ξ h iC −1 ,as − (ε(uh ), divAW ξ h ) = 0 ,
(σ AW − C ε(uAW h h ), ∇vh ) = 0
(5.32)
RT,⊥ for all ξ h ∈ Xh and vh ∈ Vh . Similarly, setting τ h = div−1 in (4.1), we RT ξ h ∈ Σh obtain −1 −1 RT RT (div σ RT h , ξ h ) + hσ h , divRT ξ h iC −1 ,as − (ε(uh ), divRT ξ h ) = 0 , RT (σ RT h − C ε(uh ), ∇vh ) = 0
(5.33)
15
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
for all ξ h ∈ Xh and vh ∈ Vh . Combining (5.32) and (5.33) leads to −1 −1 AW AW (div (σ AW − σ RT − σ RT − uRT h h ), ξ h ) + hσ h h , divAW ξ h iC −1 ,as − (ε(uh h ), divAW ξ h ) −1 −1 −1 −1 RT = hσ RT h , divRT ξ h − divAW ξ h iC −1 ,as − (ε(uh ), divRT ξ h − divAW ξ h ) , AW (σ AW − σ RT − uRT h h − C ε(uh h ), ∇vh ) = 0
or, equivalently, −1 −1 AW (div (σ AW − σ RT (PAW (σ AW − σ RT − uRT h h ), ξ h ) + (C h h )) − ∇(uh h ), divAW ξ h ) −1 −1 −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ), divRT ξ h − divAW ξ h ) + (as(C RT ξ h − divAW ξ h ) , AW −(div (σ AW − σ RT − uRT h h ), vh ) − C ε(uh h ), ε(vh )) = 0 RT for all ξ h ∈ Xh and vh ∈ Vh . In terms of the differences η h = divσ AW h −divσ h ∈ Xh AW RT and eh = uh − uh ∈ Vh , this may be written as −1 (η h , ξ h ) + (C −1 (div−1 AW η h ), divAW ξ h ) + (eh , ξ h ) −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT ξ h − divAW ξ h ) ,
(5.34)
(η h , vh ) + (C ε(eh ), ε(vh )) = 0 . Setting ξ h = η h ∈ Xh and vh = eh ∈ Vh in (5.34) leads to −1 (η h , η h ) + (C −1 (div−1 AW η h ), divAW η h ) + (eh , ξ h ) −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) ,
(η h , eh ) + (C ε(eh ), ε(eh )) = 0 . which may be combined into one equation as 2 1/2 kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k + 2(eh , η h ) + kC −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) .
(5.35)
For the left-hand side in (5.35), −1 −1 (eh , η h ) = (eh , div div−1 AW η h ) = −(∇eh , divAW η h ) = −(ε(eh ), divAW η h )
implies 2 1/2 kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k + 2(eh , η h ) + kC −1 2 1/2 = kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k − 2(ε(eh ), divAW η h ) + kC 1/2 = kη h k2 + kC −1/2 (div−1 ε(eh )k2 = G(div−1 AW η h ) − C AW η h , eh )
with the functional G defined in (2.3). For the right-hand side in (5.35), we obtain −1 RT −1 RT (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) −1 RT −1 RT ≤ kC −1 σ RT σ h )k kdiv−1 h − ε(uh )k + kas(C RT η h − divAW η h k RT 2 −1 RT 2 1/2 . h kC −1 σ RT σ h )k kη h k h − ε(uh )k + kas(C RT 1/2 . hF(σ RT kη h k , h , uh )
16
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
where the estimate −1 −1 −1 kdiv−1 η h − div−1 η h − div−1 RT η h − divAW η h k ≤ kdiv AW η h k + kdiv RT η h k . h kη h k
is used which follows again from Lemma 5.2. Therefore, (5.35) implies RT RT 1/2 G(div−1 kη h k , AW η h , eh ) . hF(σ h , uh )
which, combined with (3.9), leads to RT 1/2 kη h k . hF(σ RT h , uh )
(5.36)
and from there to (5.30). The implication of Theorems 5.3 and 5.4 is that the momentum balance error of σ RT h is of higher order compared to the error measured in the least squares functional, AW AW 1/2 RT 1/2 kdiv σ RT + F(σ RT . (5.37) h k . h F(σ h , uh ) h , uh ) This property relies on the close connection of our modified least squares formulation using Raviart-Thomas elements to the classical mixed formulation using ArnoldWinther elements. The regularity assumptions (5.15) and (5.19) are satisfied if all interior angles at points where traction boundary conditions meet displacement conditions are sufficiently small. If the problem is not H 2 regular in this sense, then, an estimate similar to (5.37) holds with h replaced by hρ for some ρ < 1. This can be shown in the usual way using interpolation arguments in scales of Sobolev spaces. The best location for the observation of (5.37) in our numerical results is in Table 4.6 where a more regular Cook’s membrane problem is treated. While the functional is reduced at a rate somewhat slower than proportional to h, the square of the momentum balance error almost behaves proportionally to h2 . We end this section with a comparison of our approach to mixed methods of saddle point type in terms of the overall dimension of the systems. We concentrate on the case k = 1 also used in our computations and denote by NP , NE , and NT the number of points, edges and triangles, respectively, in our triangulation. Then, the stress space using RT1 elements involves 4NE + 4NT unknowns. For the approach studied in our paper this needs to be augmented with the 2NP + 2NE degrees of freedom associated with the conforming quadratic elements. In the P EERS1 formulation of [12] with comparable approximation properties the RT1 elements are combined with discontinuous linear elements for the displacements (6NT unknowns) and discontinuous quadratic elements for the antisymmetry (6NT unknowns). Using the fact that for sufficiently fine triangulations NE ≈ 3NP and NT ≈ 2NP holds, we have approximately 28NP degrees of freedom in our formulation compared to 44NP in the P EERS1 approach. We thank the anonymous referees for their careful reading and valuable suggestions. REFERENCES [1] D. N. Arnold, G. Awanou, and R. Winther, Finite elements for symmetric tensors in three dimensions, Math. Comp., 77 (2008), pp. 1229–1251. [2] D. N. Arnold, F. Brezzi, and J. Douglas, PEERS: A new mixed finite element for plane elasticity, Japan J. Appl. Math., 1 (1984), pp. 347–367.
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
17
[3] D. N. Arnold and R. Winther, Mixed finite elements for elasticity, Numer. Math., 92 (2002), pp. 401–419. [4] M. Berndt, T. A. Manteuffel, and S. F. McCormick, Local error estimates and adaptive refinement for first-order system least squares, Electr. Trans. Numer. Anal., 6 (1997), pp. 35–43. [5] P. Bochev and M. Gunzburger, Least-Squares Finite Element Methods, Springer, Berlin, 2009. [6] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, Cambridge, 2nd ed., 2001. [7] J. Brandts, Y. Chen, and J. Yang, A note on least-squares mixed finite elements in relation to standard and mixed finite elements, IMA J. Numer. Anal., 26 (2006), pp. 779–789. [8] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, New York, 1991. [9] Z. Cai, J. Korsawe, and G. Starke, An adaptive least squares mixed finite element method for the stress-displacement formulation of linear elasticity, Numer. Methods Partial Differential Equations, 21 (2005), pp. 132–148. [10] Z. Cai and G. Starke, First-order system least squares for the stress-displacement formulation: Linear elasticity, SIAM J. Numer. Anal., 41 (2003), pp. 715–730. , Least squares methods for linear elasticity, SIAM J. Numer. Anal., 42 (2004), pp. 826– [11] 842. ¨ rth, On the stability of BDMS and PEERS elements, Numer. [12] M. Lonsing and R. Verfu Math., 99 (2004), pp. 131–140. [13] , A posteriori error estimators for mixed finite element methods in linear elasticity, Numer. Math., 97 (2004), pp. 757–778. [14] J. E. Pasciak and Y. Wang, A multigrid preconditioner for the mixed formulation of linear plane elasticity, SIAM J. Numer. Anal., 44 (2006), pp. 478–493. ¨ der, and G. Starke, A modified least-squares mixed finite element with [15] A. Schwarz, J. Schro improved momentum balance, Int. J. Numer. Meth. Engrg., 81 (2010), pp. 286–306. [16] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math., 53 (1988), pp. 513–538.
∗,
† ¨ ¨ ALEXANDER SCHWARZ , AND JORG SCHRODER
Abstract. A modified first-order system least squares formulation for linear elasticity, obtained by adding the antisymmetric displacement gradient in the test space, is analyzed. This approach leads to surprisingly small momentum balance error compared to standard least squares approaches. It is shown that the modified least squares formulation is well-posed and its performance is illustrated by adaptive finite element computation based on using a closely related least squares functional as a posteriori error estimator. The results of our numerical computations show that, for the modified least squares approach, the momentum balance error converges at a much faster rate than the overall error. We prove that this is due to a strong connection of the stress approximation to that obtained from a mixed formulation based on the Hellinger-Reissner principle (with exact local momentum balance). The practical significance is that our proposed approach is almost momentum-conservative at a smaller number of degrees of freedom than mixed approximations with exact local momentum balance. Key words. First-order system least squares, linear elasticity, incompressible, momentum balance, a posteriori error estimator, Raviart-Thomas elements. AMS subject classifications. 65M60, 65M15
1. Introduction. In this paper, a modified first-order system least squares formulation for linear elasticity with improved momentum balance is analyzed. Its construction is obtained from the least squares approach in [11] by adding the antisymmetric displacement gradient in the test space. This formulation was introduced in [15] from an engineering perspective and shown to give much better results in bending dominated situations. The resulting variational formulation is based on a nonsymmetric bilinear form. The ellipticity properties of this bilinear form are shown to be essentially the same as for the symmetric bilinear form associated with the “true” least squares formulation from [11]. The numerical results in [15] show a significantly improved momentum balance for our nonsymmetric formulation compared to the standard least squares method. This seems to give an explanation for the improved results, especially for bending dominated problems. In this paper, we are able to show that the reduction rate for the momentum balance error is indeed of higher order than the overall approximation error. The improvement of the reduction rate depends on the regularity of the problem. A major step in our analysis is the derivation of a close connection between the stress approximation of our method and that obtained by a mixed formulation using the Arnold-Winther element ([3], generalized to three dimensions in [1]). Such a connection was established in [7] for the least squares formulation of second-order elliptic boundary value problem and its related saddle point and primal formulations. Least squares finite element methods became very popular in recent years in many different application areas, see [5]. An important property of mixed formulations based on least squares functionals is the fact that its local evaluation provides a local a posteriori error estimate (see e.g. [4] or, in the context of linear elasticity ∗ Institut f¨ ur Angewandte Mathematik, Leibniz Universit¨ at Hannover, Welfengarten 1, 30167 Hannover, Germany ([email protected]). † Institut f¨ ur Mechanik, Fakult¨ at f¨ ur Ingenieurwissenschaften, Universit¨ at Duisburg-Essen, Universit¨ atsstr. 15, 45117 Essen, Germany ({alexander.schwarz,j.schroeder}@uni-due.de).
1
2
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
[9]). Most importantly, such approaches simultaneously construct approximations to all variables of interest (like displacements and stresses) in appropriate spaces. Alternatively, mixed methods based on saddle point formulations may be used for this purpose. For this approach, however, suitable combinations of finite element spaces are restricted by the inf-sup conditions. Families of finite element spaces for such mixed formulations of saddle point type may be found e.g. in [12]. Those families of finite element spaces are based on Raviart-Thomas elements (as a generalization of PEERS, see [2]) or Brezzi-Douglas-Marini elements (as a generalization of BDMS, see [16]). A posteriori error estimators were constructed in [13] for these saddle point approaches based on the Helmholtz decomposition. Of course, the Arnold-Winther elements from [3] which are used in this contribution as a theoretical tool could in principle also be chosen. However, these finite element spaces involve an even higher number of unknowns, in general. The outline of this paper is as follows. In the following section our modified least squares formulation for the linear elasticity model is introduced. Section 3 contains the ellipticity result for our new variational approach. We proceed with the finite element formulation and the results of our numerical computations on uniformly and adaptively refined triangulations in section 4. Finally, the strong connection between our new formulation and the mixed formulation of saddle point type based on the Arnold-Winther elements as well as the standard displacement formulation is established in section 5. 2. The Modified First-Order System Least Squares Formulation. We start from the equations of linear elasticity in the form div σ = 0 , σ − Cε(u) = 0 ,
(2.1)
where ε(u) = (∇u + ∇uT )/2 denotes the linear strain tensor and C describes the linear material law given by ν E (tr ε(u)) I ε(u) + Cε(u) = 1+ν 1 − 2ν with the elasticity modulus E and Poisson ratio ν. The boundary of Ω ⊂ IRd is divided into two parts as ∂Ω = ΓD ∪ ΓN . We consider boundary conditions of the form u = 0 on ΓD , σ · n = g on ΓN and define 1 HD (Ω) = {q ∈ H 1 (Ω) : q = 0 on ΓD }
HN (div, Ω) = {w ∈ H(div, Ω) : w · n = 0 on ΓN } . Throughout this paper, the norm on L2 (Ω) (or L2 (Ω)d , L2 (Ω)d×d , respectively) will be abbreviated by k · k and the corresponding inner product by (·, ·). We assume that ΓD is a subset of ∂Ω with positive measure (length if d = 2, area in the three-dimensional case) such that Korn’s inequality is valid in the form 1 kvk2 + k∇vk2 ≤ CK kC 1/2 ε(v)k2 for all v ∈ HD (Ω)d
(2.2)
with a constant CK > 0 (cf. [6, Section VI.3]). Note that CK is independent of ν due to the fact that kC 1/2 ε(v)k2 =
E ν E kε(u)k2 + ktr ε(u))k2 ≥ kε(u)k2 . 1+ν 1 − 2ν 1+ν
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
3
We will frequently use the decomposition of an arbitrary matrix-valued function τ ∈ L2 (Ω)d×d into its symmetric and antisymmetric part, τ = sy τ + as τ with sy τ =
τ + τT τ − τT , as τ = . 2 2
Obviously, (sy τ , as τ ) = 0 which implies kτ k2 = ksy τ k2 + kas τ k2 . If σ N ∈ H(div, Ω)d is chosen such that σ N · n = g on ΓN , then the solution of (2.1) may be obtained from minimizing the least squares functional G(σ, u) = kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2
(2.3)
1 among all σ ∈ σ N + HN (div, Ω)d and u ∈ HD (Ω)d . This is exactly the least squares functional studied in [9]. One may also extend this functional by the redundant equation
as (C −1/2 σ) = 0 (cf. [10]) which leads to F(σ, u) = kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2 + kas (C −1/2 σ)k2 .
(2.4)
The minimization of (2.4) is equivalent to the variational problem of finding σ ∈ 1 σ N + HN (div, Ω)d and u ∈ HD (Ω)d such that (div σ, div τ ) + (C −1 σ − ε(u), τ ) + (as (C −1 σ), as τ ) = 0 , (σ − C ε(u), ε(v)) = 0
(2.5)
1 (Ω)d . The additional term associated with holds for all τ ∈ HN (div, Ω)d and v ∈ HD the antisymmetric stress part is not necessary in the above least squares formulation but will be required in the following modification. Replacing ε(v) by ∇v in the displacement variation of (2.5) leads to our modified least squares method. It consists 1 in finding σ ∈ σ N + HN (div, Ω)d and u ∈ HD (Ω)d such that
(div σ, div τ ) + (C −1 σ − ε(u), τ ) + (as (C −1 σ), as τ ) = 0 , (σ − Cε(u), ∇v) = 0
(2.6)
1 holds for all τ ∈ HN (div, Ω)d and v ∈ HD (Ω)d .
3. Ellipticity. We consider the nonsymmetric bilinear form B(σ, u; τ , v) = (div σ, div τ ) + (C −1 σ − ε(u), τ − C ∇v) + (as (C −1 σ), as τ ) (3.1) associated with the modified least squares formulation (2.6). We will show that this bilinear form is coercive and continuous with respect to the product space HN (div, Ω)d × 1 HD (Ω)d equipped with the norms 1/2 kτ kdiv,C −1 = kdiv τ k2 + kC −1/2 τ k2 , 1/2 kvk1,C = kC 1/2 ∇vk2 + kvk2 .
(3.2)
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
4
For these estimates we rely on the assumption that E is on the order of one. Note that this can always be achieved by a suitable rescaling of the variables. Coercivity and continuity of the bilinear form is then achieved uniformly in the incompressible limit, i.e. as ν → 1/2. This is the statement of the following theorem. Theorem 3.1. There exist positive constants α and β which are independent of ν such that B(τ , v; τ , v) ≥ α kτ k2div,C −1 + kvk21,C (3.3) 1 holds for all (τ , v) ∈ HN (div, Ω)d × HD (Ω)d and
1/2 1/2 kτ k2div,C −1 + kvk21,C B(σ, u; τ , v) ≤ β kσk2div,C −1 + kuk21,C
(3.4)
1 holds for all (τ , v), (σ, u) ∈ HN (div, Ω)d × HD (Ω)d . Proof. The proof of (3.4) is rather straightforward. The Cauchy-Schwarz inequality immediately implies 1/2 B(σ, u; τ , v) ≤ kdiv σk2 + kC −1/2 σ − C 1/2 ε(u)k2 + kas (C −1/2 σ)k2 (3.5) 1/2 kdiv τ k2 + kC −1/2 τ − C 1/2 ∇vk2 + kas (C −1/2 τ )k2 .
Using the orthogonal splitting of an arbitrary matrix into its symmetric and antisymmetric part, we obtain kC 1/2 ∇vk2 = (C ∇v, ∇v) = (sy (C ∇v), ε(v)) + (as (C ∇v), as ∇v) = (C ε(v), ε(v)) + (C as ∇v, as ∇v) = kC 1/2 ε(v)k2 + kC 1/2 as ∇vk2 , which holds due to E ν sy ∇v + (tr ∇v) I 1+ν 1 − 2ν E ν = ε(v) + (tr ε(v)) I = Cε(v) , 1+ν 1 − 2ν E as (C ∇v) = as ∇v = C as ∇v . 1+ν
sy (C ∇v) =
In an analogous way, one obtains kC −1/2 τ k2 = kas (C −1/2 τ )k2 + ksy (C −1/2 τ )k2 . These identities imply kC 1/2 ε(v)k ≤ kC 1/2 ∇vk , kas C −1/2 τ k ≤ kC −1/2 τ k ,
(3.6)
respectively. Using (3.6), (3.5) can be further estimated as 1/2 B(σ, u; τ , v) ≤ kdiv σk2 + 2kC −1/2 σk2 + 2kC 1/2 ∇uk2 + kC −1/2 σk2 1/2 kdivτ k2 + 2kC −1/2 τ k2 + 2kC 1/2 ∇vk2 + kC −1/2 τ k2 .
(3.7)
5
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
This means that (3.4) holds with β = 3. Naturally, the coercivity estimate (3.3) is harder to establish. We start from B(τ , v; τ , v) = kdiv τ k2 + (C −1 τ − ε(v), τ − C ∇v) + kas (C −1/2 τ )k2 and note that all three terms on the right hand side are nonnegative. This is trivial for the first and third term. The second term constitutes a quadratic functional which 1 clearly attains its minimum for (τ ∗ , v∗ ) ∈ HN (div, Ω)d × HD (Ω)d satisfying 2(C −1 τ ∗ , τ ) − (ε(v∗ ) + ∇v∗ , τ ) = 0 for all τ ∈ HN (div, Ω)d , 1 2(C ε(v∗ ), ε(v)) − (τ ∗ , ε(v) + ∇v) = 0 for all v ∈ HD (Ω)d .
In particular, 2(C −1 τ ∗ , τ ∗ ) − (ε(v∗ ) + ∇v∗ , τ ∗ ) = 0 , 2(C ε(v∗ ), ε(v∗ )) − (τ ∗ , ε(v∗ ) + ∇v∗ ) = 0 which implies (C −1 τ ∗ − ε(v∗ ), τ ∗ − C ∇v∗ ) = (C −1 τ ∗ , τ ∗ ) − (τ ∗ , ∇v∗ + ε(v∗ )) + (C ε(v∗ ), v∗ ) = 0 . With a constant C1 ≥ 1 which is still free to be chosen appropriately, we therefore have 1 kdiv τ k2 + (C −1 τ − ε(v), τ − C ∇v) + C1 kas (C −1/2 τ )k2 B(τ , v; τ , v) ≥ C1 1 = kdiv τ k2 + (C −1 τ − ε(v), τ − C ε(v)) C1 (3.8) −(C −1 τ − ε(v), as (C ∇v)) + C1 kas C −1/2 τ k2 1 = kdiv τ k2 + kC −1/2 τ − C 1/2 ε(v)k2 C1 −(as (C −1/2 τ ), as (C 1/2 ∇v)) + C1 kas (C −1/2 τ )k2 . We will rely on the coercivity proof for the symmetric bilinear form from [9, Theorem 2.1] which states that, with a positive constant C0 (independently of ν), G(τ , v) ≥ C0 kdiv τ k2 + kC −1/2 τ k2 + kC 1/2 ε(v)k2 (3.9) 1 holds for all (τ , v) ∈ HN (div, Ω)d × HD (Ω)d . From Korn’s inequality (2.2) we deduce
E E kas ∇vk2 ≤ k∇vk2 1+ν 1+ν CK E CK E ν ≤ kε(v)k2 ≤ kε(v)k2 + kdiv vk2 1+ν 1+ν 1 − 2ν
kas (C 1/2 ∇v)k2 = (as (C ∇v), as ∇v) =
(3.10)
= CK (C ε(v), ε(v)) = CK kC 1/2 ε(v)k2 . Combining (3.10) with (3.9) leads to kas (C 1/2 ∇v)k2 ≤
CK G(τ , v) . C0
(3.11)
6
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
Therefore, (3.8) can be further bounded from below as follows: 1 B(τ , v; τ , v) ≥ kdivτ k2 + kC −1/2 τ − C 1/2 ε(v)k2 + C1 kas C −1/2 τ k2 C1 −(as (C −1/2 τ ), as (C 1/2 ∇v)) 1 kdiv τ k2 + kC −1/2 τ − C 1/2 ε(v)k2 + C1 kas (C −1/2 τ )k2 ≥ C1 1 1/2 2 −1/2 2 kas (C ∇v)k −C1 kas (C τ )k − 4C1 1 1 1/2 2 = G(τ , v) − kas (C ∇v)k C1 4C1 1 CK ≥ 1− G(τ , v) . C1 4C0 C1
(3.12)
If we choose C1 = max{1, CK /(2C0 )}, we have CK /(4C0 C1 ) ≤ 1/2 and therefore B(τ , v; τ , v) ≥
1 G(τ , v) . 2C1
Combined with (3.9) and (2.2), this leads to C0 1 2 −1/2 2 2 1/2 2 B(τ , v; τ , v) ≥ kdiv τ k + kC τk + kvk + kC ∇vk 2C1 CK
(3.13)
(3.14)
which finishes the proof. The constants α and β in Theorem 1 are independent of ν assuming that E is on the order of 1. This means that our coercivity and continuity estimates hold uniformly in the incompressible limit as ν → 1/2. Coercivity and continuity could also be established with respect to unscaled norms 1/2 k · kdiv := k · kdiv,I = kdiv ( · ) k2 + k · k2 , 1/2 k · k1 := k · k1,I = k∇ ( · ) k2 + k · k2 . However, these estimates would not be uniform in the incompressible limit. The fact that (3.3) and (3.4) hold ensures that our modified least squares formulation is well-posed, i.e. the variational problem (2.6) has a unique solution which depends in a continuous way on the boundary data. 4. Finite Element Approximation. For the finite element approximation of the variational problem (2.6), we choose finite-dimensional subspaces Σh ⊂ Σ := 1 HN (div, Ω)d and Vh ⊂ V := HD (Ω)d . The discrete variational problem consists in finding σ h = σ N + Σh and uh ∈ Vh such that (div σ h , div τ h ) + (C −1 σ h − ε(uh ), τ h − C ∇vh ) + (as (C −1 σ h ), as τ h ) = 0
(4.1)
holds for all τ h ∈ Σh and vh ∈ Vh . Combined with (2.6), this implies that (div (σ − σ h ), div τ h ) + (C −1 (σ − σ h ) − ε(u − uh ), τ h − C ∇vh ) + (as (C −1 (σ − σ h )), as τ h ) = 0 is satisfied for all τ h ∈ Σh and vh ∈ Vh .
(4.2)
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
7
In our computational results, next-to-lowest order Raviart-Thomas elements (RT1 ) for the stress space Σh are combined with conforming quadratic elements (P2 ) for the displacement space Vh . Example 1. We consider Cook’s membrane problem with corners (0, 0), (48, 44), (48, 60) and (0, 44) which is clamped at the left boundary segment (see figure 4.1). The rest of the boundary is subjected to surface traction forces which are zero on the upper and lower boundary segments and pointing upwards on the right boundary segment. We consider an almost incompressible material with Poisson ratio ν = 0.49.
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 8 32 128 512 2048 8192
dim Σh 40 144 544 2112 8320 33024
dim Vh 72 304 1248 5056 20352 81664
F(σ h , uh ) 12.729 5.402 2.642 1.206 0.617 0.421
kdiv σ h k2 2.640 2.125 1.240 0.566 0.387 0.334
Table 4.1 Symmetric least squares formulation (2.5): RT1 × P2 , uniform refinement
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 8 32 128 512 2048 8192
dim Σh 40 144 544 2112 8320 33024
dim Vh 72 304 1248 5056 20352 81664
F(σ h , uh ) 19.145 9.421 4.468 1.631 0.753 0.527
kdiv σ h k2 4.527 2.336 0.901 0.499 0.374 0.289
Table 4.2 Nonsymmetric least squares formulation (2.6): RT1 × P2 , uniform refinement
A comparison of the results in tables 4.1 and 4.2 shows that, on the same triangulations resulting from uniform refinement, the nonsymmetric formulation produces a faster decrease of the momentum balance error. Note that the overall least squares functional must be smaller for the symmetric formulation (on identical triangulations) since it is actually minimized. However, the higher reduction rate of the momentum balance error is not very pronounced in this example. We also compare these two formulations on a sequence of adaptively refined triangulations based on the a posteriori error estimator provided by the local evaluation of the functional F(σ h , uh ). A fixed proportion of elements with largest local contribution to the functional is refined in our adaptive approach (about 15 % of elements in our computations). The results in tables 4.3 and 4.4 show that on adaptively refined triangulations, the difference between the two formulations is much more striking. The reduction rate of the momentum balance term in the functional is significantly faster in the nonsymmetric formulation. Surprisingly, even the least squares functional on the finest refinement level is smaller for the nonsymmetric formulation. This means that these refined triangulations are better suited to our problem. The adaptively refined triangulations as well as the deformed domain are shown in figure 4.1. The most striking difference is that the nonsymmetric formulation (on
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
8
l=0 l=1 l=2 l=3 l=4 l=5 l=6 l=7 l=8 l=9 l = 10
# elem. 8 19 43 90 184 372 732 1359 2440 4351 7568
dim Σh 40 90 194 396 796 1572 3044 5586 9960 17634 30600
dim Vh 72 176 408 864 1780 3636 7204 13440 24200 43280 75352
F(σ h , uh ) 12.729 5.476 2.584 1.627 1.061 0.748 0.519 0.413 0.341 0.267 0.185
kdiv σ h k2 2.640 2.048 1.257 0.975 0.718 0.543 0.394 0.323 0.253 0.162 0.082
Table 4.3 Symmetric least squares formulation (2.5): RT1 × P2 , adaptive refinement
l=0 l=1 l=2 l=3 l=4 l=5 l=6 l=7 l=8 l=9 l = 10
# elem. 8 19 41 87 190 370 710 1346 2464 4415 7774
dim Σh 40 90 186 382 812 1560 2948 5528 10044 17902 31408
dim Vh 72 176 388 836 1848 3620 6992 13316 24452 43908 77428
F(σ h , uh ) 19.145 10.773 4.595 2.558 1.277 0.796 0.592 0.476 0.360 0.230 0.129
kdiv σ h k2 4.527 3.361 1.874 1.010 0.548 0.379 0.269 0.151 0.065 0.020 0.005
Table 4.4 Nonsymmetric least squares formulation (2.6): RT1 × P2 , adaptive refinement
the right) puts much more emphasis at refining near the upper left corner. Example 2. We consider the more regular Cook’s membrane problem with corners (0, 0), (48, 14), (48, 30) and (0, 44) where the singularities are much less severe. This is a good illustration of our analysis in the next section which will rely on the regularity of the underlying boundary value problem. It will indicate that the decrease of the momentum balance error relative to the overall error will be more pronounced the more regular the problem is. We repeat our computations for this more regular modification of Cook’s membrane on a sequence of uniformly refined triangulations. The results in table 4.6 for the nonsymmetric formulation clearly show a faster reduction of the momentum balance error which the numbers in table 4.5 for the symmetric formulation do not. 5. A Strong Connection To Mixed and Displacement Formulations. In [7], Brandts, Chen and Yang have established a strong connection for second-order elliptic problems between the least squares formulation and corresponding mixed and standard Galerkin approaches. Such a strong connection also holds for the linear elasticity problem between our modified variational formulation (2.6) and an appropriate mixed and the standard displacement-based approach. This result is the basis for
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
9
Fig. 4.1. Triangulation of deformed Cook’s membrane after seven adaptive refinement steps: symmetric (left) vs. nonsymmetric (right) formulation
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 11 44 176 704 2816 11264
dim Σh 54 196 744 2896 11424 45376
dim Vh 100 420 1720 6960 28000 112320
F(σ h , uh ) 4.575 1.890 0.892 0.511 0.338 0.204
kdiv σ h k2 1.491 0.789 0.467 0.330 0.197 0.079
Table 4.5 Symmetric least squares formulation (2.5): RT1 × P2 , uniform refinement
l=0 l=1 l=2 l=3 l=4 l=5
# elem. 11 44 176 704 2816 11264
dim Σh 54 196 744 2896 11424 45376
dim Vh 100 420 1720 6960 28000 112320
F(σ h , uh ) 7.662 2.551 1.151 0.724 0.487 0.264
kdiv σ h k2 3.149 0.916 0.397 0.180 0.056 0.011
Table 4.6 Nonsymmetric least squares formulation (2.6): RT1 × P2 , uniform refinement
our proof of the improved momentum balance observed in the computational results shown in Section 4. In our analysis, we will use the finite element spaces introduced by Arnold and
10
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
Winther in [3] which are defined as follows for k ≥ 1: 2×2 ΣAW = {τ ∈ Σ : as τ = 0 , τ |T ∈ Pk+2 , div τ |T ∈ Pk2 } , h
(5.1)
where Pk denotes the space of polynomials of degree k (in two dimensions). Important for our consideration is the fact that div ΣAW coincides with the space of piecewise h polynomials of degree k (without any inter-element continuity conditions). We restrict our analysis to two dimensional problems. For some fixed integer k ≥ 1, let ΣAW ⊆ Σ denote the Arnold-Winther space of h degree k introduced above. Moreover, let Vh ⊂ V be the standard H 1 (Ω)-conforming finite element space of piecewise polynomials of degree k + 1. Finally, denote by Xh ⊂ L2 (Ω)2 the piecewise polynomial spaces of degree k without any continuity requirements. Let us denote the finite element approximation based on the variational formuAW AW AW N lation (4.1) with Σh = ΣAW by (σ AW + h h , uh ). Note that (σ h , uh ) ∈ (σ AW Σh ) × Vh satisfies as σ h = 0 and therefore also minimizes the least squares functional G(σ h , uh ) with respect to ΣAW h ) × Vh . The first step of our analysis consists AW in establishing a close connection of (σ AW h , uh ) to the following mixed and standard displacement approaches. The mixed formulation that we consider is based on Hellinger-Reissner principle and uses the finite element space by Arnold and Winther [3]. The mixed finite element m AW N approximation (σ m h , ξ h ) ∈ (σ + Σh ) × Xh with these spaces satisfy m (C −1 (σ − σ m h ), τ h ) + (u − ξ h , div τ h ) = 0
(5.2)
for all τ h ∈ ΣAW and h (div(σ − σ m h ), η h ) = 0 for all η h ∈ Xh .
(5.3)
For the combination ΣAW × Xh introduced above (Arnold-Winther space combined h with piecewise polynomials of degree k) the first constraint in (5.3) implies that 2 div σ m h = 0 holds exactly in L (Ω). Finally, the standard displacement formulation udh ∈ Vh satisfies (Cε(u − udh ), ε(vh )) = 0
(5.4)
for all vh ∈ Vh . We are now ready to prove the generalization of the Brandts-Chen-Yang formula (cf. [7, Lemma 3.2]) to our modified least squares formulation (2.6). Lemma 5.1. For some integer k ≥ 1, let ΣAW denote the Arnold-Winther space h of degree k and let Vh denote the standard continuous elements of degree k + 1. For AW AW N the finite element approximation (σ AW h , uh ) ∈ (σ + Σh ) × Vh based on (4.1), the AW m N mixed approximation σ h ∈ σ +Σh based on (5.2) and (5.3), and the displacement approximation udh ∈ Vh based on (5.4), the following relation holds: AW AW B(σ AW − σm − udh ; σ AW − σm − udh ) h h , uh h h , uh AW d AW = (C −1 (σ − σ m − σm − σm h ), σ h h ) − (ε(u − uh ), σ h h ).
(5.5)
11
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
Proof. The definition of the bilinear form (3.1) leads to AW AW B(σ AW − σm − udh ; σ AW − σm − udh ) h h , uh h h , uh 2 = kdiv (σ AW − σm h h )k AW AW + (C −1 (σ AW − σm − udh ), σ AW − σm − udh )) h h ) − ε(uh h h − C ∇(uh AW + (as (σ AW − σm − σm h h ), as (σ h h ))
= (div (σ AW − σ), div (σ AW − σm h h h ) AW + (C −1 (σ AW − σ) − ε(uAW − u), σ AW − σm − udh )) h h h h − C ∇(uh
(5.6)
+ (as (σ AW − σ), as (σ AW − σm h h h )) d AW AW + (C −1 (σ − σ m − σm − udh )) h ) − ε(u − uh ), σ h h − C ∇(uh AW + (as (σ − σ m − σm h ), as (σ h h )) AW m AW = (C −1 (σ − σ m − σm − udh )) h ), σ h h ) − (σ − σ h , ∇(uh d AW − (ε(u − udh ), σ AW − σm − udh )) h h ) + (C ε(u − uh ), ε(uh AW + (as (σ − σ m − σm h ), as (σ h h )) ,
AW m in the last equality (note that σ AW where we used (4.2) with Σh = ΣAW h h −σ h ∈ Σh and uAW − udh ∈ Vh ). The second and fourth terms in (5.6) vanish due to (5.3) and h (5.4), respectively. The last term vanishes since symmetry is imposed in the stress space ΣAW h . This leads us to the desired generalization of the Brandts-Chen-Yang formula (5.5). In the sequel, we will use discrete inverses of the divergence operator with respect to the spaces ΣRT ⊆ Σ and ΣAW ⊆ Σ (cf. [14, Sect. 2]). To this end, for any h h AW to be the solution of the saddle point problem ∈ Σ ξ h ∈ Xh , we may define ΞAW h h
(C −1 ΞAW h , τ h ) − (ξ h , div τ h ) = 0
for all τ h ∈ ΣAW , h
(div ΞAW h , η h ) = (ξ h , η h )
for all η h ∈ Xh .
(5.7)
RT Similarly, we may define ΞRT h ∈ Σh to be the solution of −1 RT Ξh ), as τ h ) − (ξ h , div τ h ) = 0 for all τ h ∈ ΣRT (C −1 ΞRT h , h , τ h ) + (as(C
(div ΞRT h , η h ) = (ξ h , η h ) for all η h ∈ Xh .
(5.8)
If we define subspaces RT/AW,0
Σh
RT/AW
= {τ h ∈ Σh
: div τ h = 0}
(5.9)
and the inner product h · , · iC −1 ,as on Σ given by h σ , τ iC −1 ,as = (C −1 σ, τ ) + (as(C −1 σ), as τ ) ,
(5.10)
AW,0 we obtain hΞAW and, similarly, hΞRT h , τ h iC −1 ,as = 0 for all τ h ∈ Σh h , τ h iC −1 ,as = 0 RT,0 for all τ h ∈ Σh . Furthermore, we define orthogonal subspaces RT/AW,⊥
Σh
RT/AW
= {τ h ∈ Σh
RT/AW,0
: hτ h , ω h iC −1 ,as = 0 for all ω h ∈ Σh
}.
(5.11)
The construction in (5.7) and (5.8) provides discrete inverses of the divergence operRT/AW,⊥ ator mapping Σh , respectively, onto Xh . We therefore use the notation −1 AW ΞRT = div−1 h = divRT ξ h and Ξh AW ξ h .
(5.12)
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
12
For any ξ ∈ L2 (Ω)2 , we may proceed in any similar way as in (5.7) and (5.8) with respect to the entire Sobolev space Σ = HN (div, Ω)2 . To this end, let Ξ ∈ Σ be the solution of (C −1 Ξ, τ ) + (as(C −1 Ξ), as τ ) − (ξ, div τ ) = 0 (div Ξ, η) = (ξ, η)
for all τ ∈ Σ , for all η ∈ L2 (Ω)2 .
(5.13)
This defines the inverse of the divergence operator mapping the orthogonal subspace Σ⊥ of Σ with respect to Σ0 = {τ ∈ Σ : div τ = 0}
(5.14)
onto L2 (Ω)2 , Ξ = div−1 ξ. From now on, we will use the notation ιh . κh to denote that ιh ≤ Cκh holds with some constant C independently of h. Some fundamental properties of the operators introduced in (5.7) and (5.8) are summarized in the following lemma. The statements of Lemma 5.2 are well-known and, in the case of the Arnold-Winther spaces, also used in [14]. We state them explicitly here for the sake of better readability. Lemma 5.2. Under the assumptions that the problem (5.13) satisfies the regularity assumption Ξ ∈ H 1 (Ω)2×2 with kΞk1,Ω . kξk ,
(5.15)
the following holds for the discrete inverses of the divergence operator defined in (5.7) and (5.8): −1 kdiv−1 ξ h − div−1 RT ξ h k . h kξ h k , kdivRT ξ h k . kξ h k
(5.16)
for all ξ h ∈ Xh and −1 kdiv−1 ξ h − div−1 AW ξ h k . h kξ h k , kdivAW ξ h k . kξ h k
(5.17)
for all ξ h ∈ Xh . Proof. The regularity assumption (5.15) implies kΞk1,Ω . kdiv Ξk .
(5.18)
The standard approximation properties of Raviart-Thomas spaces (see [8, Prop. 3.9]) lead to kΞ − ΞRT h k . h kΞk1,Ω which, combined with (5.18), proves (5.16). Finally, the first approximation property stated in [3, Theorem 6.1] implies kΞ − ΞAW h k . h kΞk1,Ω which, combined with (5.18), leads to (5.17). The next step consists in showing that the momentum balance error of the formulation (4.1) using Arnold-Winther elements converges at a faster rate than the error norm measured by the least squares functional.
13
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
m AW d Theorem 5.3. Let σ AW h , σ h , uh , uh be as in Lemma 5.1 and assume that (5.15) is satisfied. Furthermore, assume that the boundary value problem associated with (2.1) is H 2 -regular in the sense that, for any g ∈ H 1/2 (ΓN )2 ,
u ∈ H 2 (Ω)2 with kuk2,Ω . kgk1/2,ΓN
(5.19)
holds. Then, 1/2 −1/2 2 1/2 kdiv σ AW (σ − σ m (ε(u) − ε(udh ))k2 h k . h kC h )k + kC
(5.20)
and AW AW 1/2 kdiv σ AW . h k . h F(σ h , uh )
(5.21)
Proof. From (5.5) and (3.3), we obtain (note that div σ m h = 0) 2 AW 2 kdiv σ AW − σm h k = kdiv (σ h h )k AW d AW . (C −1 (σ − σ m − σm − σm h ), σ h h ) − (ε(u − uh ), σ h h ).
(5.22)
For the first term in (5.22), (5.2) leads to m AW AW (C −1 (σ − σ m − σm − σm h ), σ h h ) = (u − ξ h , div (σ h h )) m AW AW = (u − ξ m h , div σ h ) = (u − ξ h , div Ξh )
= (C
−1
(σ −
AW σm h ), Ξh )
(5.23)
,
= is the solution of (5.7) with ξ h = div σ AW (note that div ΞAW ∈ ΣAW where ΞAW h h h h AW div σ h holds). We also have −1 −1 (C −1 (σ − σ m (σ − σ m (σ − σ m h ), Ξ) = (C h ), Ξ) + (as (C h )), as Ξ) = 0 ,
(5.24)
if Ξ ∈ Σ is the solution of (5.13) with ξ = div σ AW (note that div (σ − σ m h h ) = 0 and m as (σ − σ h ) = 0 holds). The combination of (5.23) and (5.24) leads to AW AW −1 − Ξ) (C −1 (σ − σ m − σm (σ − σ m h ), σ h h ) = (C h ), Ξh −1/2 ≤ kC −1/2 (σ − σ m (Ξ − ΞAW h )kkC h )k
. h kC
−1/2
(σ −
σm h )k
kdiv
σ AW h k
(5.25)
,
AW where kC −1/2 (Ξ − ΞAW h )k . hkdiv σ h k is used which holds due to Lemma 5.2. For the second term in (5.22), integration by parts gives d AW d AW −(ε(u − udh ), σ AW − σm − σm h h ) = (u − uh , div (σ h h ) = (u − uh , div σ h ) ,
where as (σ AW − σm h h ) = 0 was used. This leads to d AW −(ε(u − udh ), σ AW − σm h h ) ≤ ku − uh k kdiv σ h k
. hkC 1/2 ε(u − udh )k kdiv σ AW h k,
(5.26)
where the estimate ku−udh k . hkC 1/2 ε(u−udh )k was used which holds due to the classical Aubin-Nitsche duality argument and (2.2) under our assumption of H 2 regularity for u. The combination of (5.25) and (5.26) proves (5.20).
14
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
The estimate (5.21) is obtained from (5.20) in the following way: The relation (5.2) implies 2 −1/2 kC −1/2 (σ − σ m (σ − τ h )k2 h )k ≤ kC
for all τ h ∈ ΣAW,0 and therefore h 2 −1/2 AW 2 kC −1/2 (σ − σ m (σ − (σ AW − div−1 h )k ≤ kC h AW div σ h )k −1 2 AW 2 ≤ 2 kC −1/2 (σ − σ AW )k + kdiv div σ )k h h AW 2 AW 2 . 2 kC −1/2 (σ − σ AW h )k + kdiv σ h k 2 AW 2 , = 2 kC −1/2 (σ − σ AW )k + kdiv (σ − σ )k h h
(5.27)
where Lemma 5.2 is used again. Similarly, the relation (5.4) leads to kC 1/2 (ε(u) − ε(udh ))k2 ≤ kC 1/2 (ε(u) − ε(vh ))k2 for all vh ∈ Vh which implies 2 kC 1/2 (ε(u) − ε(udh ))k2 ≤ kC 1/2 (ε(u) − ε(uAW h ))k .
(5.28)
From (5.27), (5.28), together with (3.9), 2 1/2 AW AW AW kC −1/2 (σ − σ m (ε(u) − ε(udh ))k2 . G(σ AW h )k + kC h , uh ) = F(σ h , uh ) (5.29)
is obtained which finishes the proof of (5.21). RT AW AW N RT N Theorem 5.4. Let (σ RT h , uh ) ∈ (σ + Σh ) × Vh and (σ h , uh ) ∈ (σ + AW Σh ) × Vh be the approximations of our modified least squares method (4.1) in the respective finite element spaces. Then, under the assumption that (5.15) holds, AW RT RT 1/2 kdiv σ RT . h k . kdiv σ h k + h F(σ h , uh )
(5.30)
Proof. We start with the observation that divσ AW as well as divσ RT h h are contained in Xh . Let AW ⊥ div−1 = PAW σ AW , h AW div σ h RT ⊥ RT div−1 RT div σ h = PRT σ h
(5.31)
be the orthogonal projections with respect to h · , · iC −1 ,as onto ΣAW,⊥ and ΣRT,⊥ , h h AW,⊥ −1 respectively. Plugging the special test functions τ h = divAW ξ h ∈ Σh with ξh ∈ Xh into (4.1) leads to −1 −1 AW AW (div σ AW h , ξ h ) + hσ h , divAW ξ h iC −1 ,as − (ε(uh ), divAW ξ h ) = 0 ,
(σ AW − C ε(uAW h h ), ∇vh ) = 0
(5.32)
RT,⊥ for all ξ h ∈ Xh and vh ∈ Vh . Similarly, setting τ h = div−1 in (4.1), we RT ξ h ∈ Σh obtain −1 −1 RT RT (div σ RT h , ξ h ) + hσ h , divRT ξ h iC −1 ,as − (ε(uh ), divRT ξ h ) = 0 , RT (σ RT h − C ε(uh ), ∇vh ) = 0
(5.33)
15
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
for all ξ h ∈ Xh and vh ∈ Vh . Combining (5.32) and (5.33) leads to −1 −1 AW AW (div (σ AW − σ RT − σ RT − uRT h h ), ξ h ) + hσ h h , divAW ξ h iC −1 ,as − (ε(uh h ), divAW ξ h ) −1 −1 −1 −1 RT = hσ RT h , divRT ξ h − divAW ξ h iC −1 ,as − (ε(uh ), divRT ξ h − divAW ξ h ) , AW (σ AW − σ RT − uRT h h − C ε(uh h ), ∇vh ) = 0
or, equivalently, −1 −1 AW (div (σ AW − σ RT (PAW (σ AW − σ RT − uRT h h ), ξ h ) + (C h h )) − ∇(uh h ), divAW ξ h ) −1 −1 −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ), divRT ξ h − divAW ξ h ) + (as(C RT ξ h − divAW ξ h ) , AW −(div (σ AW − σ RT − uRT h h ), vh ) − C ε(uh h ), ε(vh )) = 0 RT for all ξ h ∈ Xh and vh ∈ Vh . In terms of the differences η h = divσ AW h −divσ h ∈ Xh AW RT and eh = uh − uh ∈ Vh , this may be written as −1 (η h , ξ h ) + (C −1 (div−1 AW η h ), divAW ξ h ) + (eh , ξ h ) −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT ξ h − divAW ξ h ) ,
(5.34)
(η h , vh ) + (C ε(eh ), ε(vh )) = 0 . Setting ξ h = η h ∈ Xh and vh = eh ∈ Vh in (5.34) leads to −1 (η h , η h ) + (C −1 (div−1 AW η h ), divAW η h ) + (eh , ξ h ) −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) ,
(η h , eh ) + (C ε(eh ), ε(eh )) = 0 . which may be combined into one equation as 2 1/2 kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k + 2(eh , η h ) + kC −1 RT −1 RT = (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) .
(5.35)
For the left-hand side in (5.35), −1 −1 (eh , η h ) = (eh , div div−1 AW η h ) = −(∇eh , divAW η h ) = −(ε(eh ), divAW η h )
implies 2 1/2 kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k + 2(eh , η h ) + kC −1 2 1/2 = kη h k2 + kC −1/2 (div−1 ε(eh )k2 AW η h )k − 2(ε(eh ), divAW η h ) + kC 1/2 = kη h k2 + kC −1/2 (div−1 ε(eh )k2 = G(div−1 AW η h ) − C AW η h , eh )
with the functional G defined in (2.3). For the right-hand side in (5.35), we obtain −1 RT −1 RT (C −1 σ RT σ h ), div−1 h − ε(uh ) + as(C RT η h − divAW η h ) −1 RT −1 RT ≤ kC −1 σ RT σ h )k kdiv−1 h − ε(uh )k + kas(C RT η h − divAW η h k RT 2 −1 RT 2 1/2 . h kC −1 σ RT σ h )k kη h k h − ε(uh )k + kas(C RT 1/2 . hF(σ RT kη h k , h , uh )
16
¨ G. STARKE, A. SCHWARZ AND J. SCHRODER
where the estimate −1 −1 −1 kdiv−1 η h − div−1 η h − div−1 RT η h − divAW η h k ≤ kdiv AW η h k + kdiv RT η h k . h kη h k
is used which follows again from Lemma 5.2. Therefore, (5.35) implies RT RT 1/2 G(div−1 kη h k , AW η h , eh ) . hF(σ h , uh )
which, combined with (3.9), leads to RT 1/2 kη h k . hF(σ RT h , uh )
(5.36)
and from there to (5.30). The implication of Theorems 5.3 and 5.4 is that the momentum balance error of σ RT h is of higher order compared to the error measured in the least squares functional, AW AW 1/2 RT 1/2 kdiv σ RT + F(σ RT . (5.37) h k . h F(σ h , uh ) h , uh ) This property relies on the close connection of our modified least squares formulation using Raviart-Thomas elements to the classical mixed formulation using ArnoldWinther elements. The regularity assumptions (5.15) and (5.19) are satisfied if all interior angles at points where traction boundary conditions meet displacement conditions are sufficiently small. If the problem is not H 2 regular in this sense, then, an estimate similar to (5.37) holds with h replaced by hρ for some ρ < 1. This can be shown in the usual way using interpolation arguments in scales of Sobolev spaces. The best location for the observation of (5.37) in our numerical results is in Table 4.6 where a more regular Cook’s membrane problem is treated. While the functional is reduced at a rate somewhat slower than proportional to h, the square of the momentum balance error almost behaves proportionally to h2 . We end this section with a comparison of our approach to mixed methods of saddle point type in terms of the overall dimension of the systems. We concentrate on the case k = 1 also used in our computations and denote by NP , NE , and NT the number of points, edges and triangles, respectively, in our triangulation. Then, the stress space using RT1 elements involves 4NE + 4NT unknowns. For the approach studied in our paper this needs to be augmented with the 2NP + 2NE degrees of freedom associated with the conforming quadratic elements. In the P EERS1 formulation of [12] with comparable approximation properties the RT1 elements are combined with discontinuous linear elements for the displacements (6NT unknowns) and discontinuous quadratic elements for the antisymmetry (6NT unknowns). Using the fact that for sufficiently fine triangulations NE ≈ 3NP and NT ≈ 2NP holds, we have approximately 28NP degrees of freedom in our formulation compared to 44NP in the P EERS1 approach. We thank the anonymous referees for their careful reading and valuable suggestions. REFERENCES [1] D. N. Arnold, G. Awanou, and R. Winther, Finite elements for symmetric tensors in three dimensions, Math. Comp., 77 (2008), pp. 1229–1251. [2] D. N. Arnold, F. Brezzi, and J. Douglas, PEERS: A new mixed finite element for plane elasticity, Japan J. Appl. Math., 1 (1984), pp. 347–367.
LEAST SQUARES WITH IMPROVED MOMENTUM BALANCE
17
[3] D. N. Arnold and R. Winther, Mixed finite elements for elasticity, Numer. Math., 92 (2002), pp. 401–419. [4] M. Berndt, T. A. Manteuffel, and S. F. McCormick, Local error estimates and adaptive refinement for first-order system least squares, Electr. Trans. Numer. Anal., 6 (1997), pp. 35–43. [5] P. Bochev and M. Gunzburger, Least-Squares Finite Element Methods, Springer, Berlin, 2009. [6] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, Cambridge, 2nd ed., 2001. [7] J. Brandts, Y. Chen, and J. Yang, A note on least-squares mixed finite elements in relation to standard and mixed finite elements, IMA J. Numer. Anal., 26 (2006), pp. 779–789. [8] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, New York, 1991. [9] Z. Cai, J. Korsawe, and G. Starke, An adaptive least squares mixed finite element method for the stress-displacement formulation of linear elasticity, Numer. Methods Partial Differential Equations, 21 (2005), pp. 132–148. [10] Z. Cai and G. Starke, First-order system least squares for the stress-displacement formulation: Linear elasticity, SIAM J. Numer. Anal., 41 (2003), pp. 715–730. , Least squares methods for linear elasticity, SIAM J. Numer. Anal., 42 (2004), pp. 826– [11] 842. ¨ rth, On the stability of BDMS and PEERS elements, Numer. [12] M. Lonsing and R. Verfu Math., 99 (2004), pp. 131–140. [13] , A posteriori error estimators for mixed finite element methods in linear elasticity, Numer. Math., 97 (2004), pp. 757–778. [14] J. E. Pasciak and Y. Wang, A multigrid preconditioner for the mixed formulation of linear plane elasticity, SIAM J. Numer. Anal., 44 (2006), pp. 478–493. ¨ der, and G. Starke, A modified least-squares mixed finite element with [15] A. Schwarz, J. Schro improved momentum balance, Int. J. Numer. Meth. Engrg., 81 (2010), pp. 286–306. [16] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math., 53 (1988), pp. 513–538.
