for the Stress-Pressure-Displacement Elasticity
Equations***
Suh–Yuh Yang1
1Department of Applied Mathematics
National Chiao Tung University Hsinchu 30050, Taiwan
Ching L. Chang2
2Department of Mathematics
Cleveland State University Cleveland, Ohio 44115
Received May 12, 1997; accepted November 4, 1997
A new stress-pressure-displacement formulation for the planar elasticity equations is proposed by introducing the auxiliary variables, stresses, and pressure. The resulting first-order system involves a nonnegative parameter that measures the material compressibility for the elastic body. A two-stage least-squares finite element procedure is introduced for approximating the solution to this system with appropriate boundary conditions. It is shown that the two-stage least-squares scheme is stable and, with respect to the order of approximation for smooth exact solutions, the rates of convergence of the approximations for all the unknowns are optimal both in the H1-norm and in the L2-norm. Numerical experiments with various
values of the parameter are examined, which demonstrate the theoretical estimates. Among other things, computational results indicate that the behavior of convergence is uniform in the nonnegative parameter.
c
1998 John Wiley & Sons, Inc. Numer Methods Partial Differential Eq 14: 297–315, 1998
Keywords: elasticity equations; least-squares; finite elements; error estimates I. INTRODUCTION
In the last ten years, the least squares finite element techniques have been extensively applied in many different fields such as fluid dynamics [1–12], elasticity [13–17], electromagnetism [18– 19], and semiconductor device physics [20] (see also [21–23] and many references therein). The
Correspondence to: Suh–Yuh Yang
∗This work was done when C. L. Chang visited the Department of Applied Mathematics of the National Chiao Tung
University, Taiwan, in the summer of 1996. c
least-squares finite element approach represents a fairly general methodology that can produce a variety of algorithms. For the elasticity problem, Franca et al. [15, 16] proposed some finite element methods that are constructed by adding various least-squares terms to the classical mixed formulation. These methods can be subdivided into two categories for attaining stability, depend-ing on whether the Babuˇska–Brezzi condition is circumvented or satisfied. Results for the full least-squares finite element methods applied to a stress-pressure-displacement formulation with the displacement boundary conditions have been recently reported in [17]. The present article investigates a two-stage least-squares finite element procedure for treating the elasticity equations with more general boundary conditions.
Introducing the auxiliary variables (stresses and pressure), we can recast the original two-dimensional elasticity system of second-order equations as an equivalent parameter-dependent first-order system with eight equations and six unknowns, in which the nonnegative parameter measures the material compressibility for the elastic body. This new stress-pressure-displacement formulation is different from that in [15, 24, 25] but is similar to the third formulation introduced in [16]. It can be further decomposed into two dependent subsystems, the stress-pressure system and the displacement system recovered from the stresses and pressure. Moreover, we can prove that the stress-pressure system with appropriate boundary conditions is an elliptic system in the sense of Petrovski and satisfies the Lopatinski condition [26]. Taking advantage of these properties, we propose a two-stage least-squares finite element procedure for these two subsystems to obtain approximations of all the unknowns in an orderly way.
The two-stage least-squares finite element presented approach offers many advantages:
•
Because the two-stage procedure leads to two minimization problems (rather than the saddle point problem resulting from the mixed finite element procedure), the approximation spaces need not satisfy the Babuˇska–Brezzi condition, and a single continuous piecewise polynomial space can be used for approximating all the unknowns in both stages.•
Its discretization results in two symmetric and positive definite linear algebraic systems both with condition number O(h−2), where h is the mesh parameter. This allows the useof efficient solvers such as the conjugate gradient method to solve the corresponding large linear systems.
•
Accurate approximations of the stresses, pressure, and the displacements can be obtained in an orderly way according to the two-stage procedure.•
Under suitable regularity assumptions, the least-squares approximations for all the un-knowns have optimal order of approximation in the H1-norm and in the L2-norm.•
Numerical experiments with various values of the parameter are examined, which confirm the theoretical error estimates. Among other things, computational results indicate that the behavior of convergence is uniform in the nonnegative parameter.In addition, compared with the least-squares finite element methods developed in [17], the most significant features of the present approach are the following:
•
The methods used in [17] work well only for the displacement boundary conditions, but the present two-stage least-squares method can be applied to the general stress-pressure-displacement boundary conditions, which are more useful in practical applications.•
Just as for the two-stage methods [5, 27], the proposed method has a computational advan-tage over the methods in [17]. Indeed, a linear system of size 6N must be solved for the methods in [17], where N is the dimension of the common approximation space. However, the two-stage procedure requires only the solution of a system of size 4N, followed by the solution of a system of size 2N, each with smaller bandwidths.The remainder of the article is organized as follows. In Section II, we propose a new stress-pressure-displacement formulation for the elasticity equations. This is then decomposed into two subsystems, the stress-pressure system and the displacement system, with respective appropriate boundary conditions. In Section III, a two-stage least-squares finite element procedure is given, as well as its fundamental properties. In Section IV, a priori estimates for the stress-pressure system are derived. In Section V, error analysis is presented. In Section VI, the condition numbers of the resulting linear systems are estimated. Finally, in Section VII, some numerical experiments are examined to demonstrate this approach.
II. PRELIMINARIES
We shall consider the numerical solution of the boundary value problem,
−2µ ∇ · ε(u) + 1 − 2νν ∇(∇ · u) = f in Ω, (2.1) u = 0 on Γ1, (2.2) 2µ ε(u) + 1 − 2νν (∇ · u)I · n = g on Γ2, (2.3)
with the following notation:
•
Ω ⊂ R2is a bounded domain representing the region occupied by an elastic body.•
Γ := ∂Ω is the smooth boundary of Ω, which is partitioned into two disjoint open parts, Γ1and Γ2, such that Γ = ¯Γ1∪ ¯Γ2and measure (Γ1) > 0.•
µ is the shear modulus given byµ = 2(1 + ν)E > 0,
where ν is the Poisson ratio, 0 < ν < 0.5, and E > 0 is the Young modulus. The upper limit of the Poisson ratio, i.e., ν → 0.5−, corresponds to an incompressible material.
•
u = (u1, u2)tis the displacement field.•
f = (f1, f2)tis the density of a body force acting on the body.•
g = (g1, g2)tis the density of a surface force acting on Γ2.•
n = (n1, n2)tis the outward unit normal vector to ∂Ω.•
ε(u) is the strain tensor given byε(u) = (εij(u))2×2= (12(∂jui+ ∂iuj))2×2.
•
I is the 2 × 2 identity matrix.Introducing the auxiliary variables, ϕ1, ϕ2, ϕ3, and p, such that
ϕ1=∂u∂x1, (2.4)
ϕ3=∂u∂x2, (2.6) −ϕ1−1 − 2νν p = ∂u∂y2 (2.7) on ¯Ω and letting = 1 − 2νν > 0, 0 < ν < 12, we can rewrite (2.1) as 2µ −∂ϕ1 ∂x − 1 2 ∂ϕ2 ∂y − 1 2 ∂ϕ3 ∂y + ∂p ∂x = f1 in Ω, (2.8) 2µ ∂ϕ1 ∂y − 1 2 ∂ϕ2 ∂x − 1 2 ∂ϕ3 ∂x + (1 + ) ∂p ∂y = f2 in Ω. (2.9)
We call ϕi the stresses and p the artificial pressure, and remark that the ‘‘pressure’’ p gives
the hydrostatic pressure only in the incompressible limit (cf. Remark 2.1 below). Note that a combination of ϕ1, ϕ2, ϕ3and p can represent the actual stresses σij(for i, j = 1, 2), which are
given by σ(u) = (σij(u))2×2= 2µ ε(u) + ν 1 − 2ν(∇ · u)I .
Also, by (2.4)–(2.7), we obtain the following two compatibility equations:
∂ϕ1 ∂y − ∂ϕ2 ∂x = 0 in Ω, (2.10) ∂ϕ1 ∂x + ∂ϕ3 ∂y + ∂p ∂x = 0 in Ω. (2.11)
System of Eqs. (2.4)–(2.11) is the so-called stress-pressure-displacement formulation for the two-dimensional elasticity equations, which is different from those in [15, 24, 25] but is similar to the third formulation introduced in [16]. Moreover, we can show that a sufficiently smooth solution of (2.1) solves system (2.4)–(2.11), and vice versa. It is interesting to observe that the relations between the stresses ϕi, pressure p, and the displacements uiare defined by
Eqs. (2.4)–(2.7), and the stress-pressure system (2.8)–(2.11) is independent of the displacements
ui. Therefore, if one can solve (2.8)–(2.11) with appropriate boundary conditions, then the displacements can be recovered from the stresses and pressure by solving Eqs. (2.4)–(2.7) with the boundary requirement u = 0 on Γ1. Our two-stage procedure is thus motivated.
Remark 2.1. For the incompressible limit, = 0, the first-order system (2.4)–(2.11) is the system of stress-pressure-velocity Stokes equations, which have been studied in [5, 8]. In the context, u represents the velocity field for the Stokes flow, p expresses the pressure with appropriate scaling, and µ denotes the inverse of the Reynolds number. We also remark that all the results developed below still hold for the case of = 0.
To deal with the boundary conditions, we note that (2.2) implies that the tangential derivatives of uivanish, ∇ui· (n2, −n1)t= 0, i = 1, 2, that is,
n1ϕ1+ n2ϕ3+ n1p = 0 on Γ1. (2.13) Combining (2.12)–(2.13) with
n1u1+ n2u2= 0 on Γ1,
we can verify (2.2) as well. Also, boundary conditions (2.3) can be written as
2µn1ϕ1+ µn2ϕ2+ µn2ϕ3− 2µn1p = g1 on Γ2. (2.14)
−2µn2ϕ1+ µn1ϕ2+ µn1ϕ3− 2µ(1 + )n2p = g2 on Γ2. (2.15) It is now clear that our strategy is to solve the stress-pressure system (2.8)–(2.11) with boundary conditions (2.12)–(2.15) at the first stage, i.e., to solve
LspΦ := AspΦx+ BspΦy= F in Ω, (2.16) RspΦ := CspΦ = G on Γ, (2.17) where Φ = (ϕ1, ϕ2, ϕ3, p)t, F = (f1, f2, 0, 0)t, Asp= −2µ 0 0 2µ 0 −µ −µ 0 0 −1 0 0 1 0 0 , Bsp= 0 −µ −µ 0 2µ 0 0 2µ(1 + ) 1 0 0 0 0 0 1 0 , Csp= n2 −n1 0 0 n1 0 n2 n1 and G = 0 0 on Γ1, Csp= 2µn1 µn2 µn2 −2µn1 −2µn2 µn1 µn1 −2µ(1 + )n2 and G = g1 g2 on Γ2.
The second stage is to solve the displacement system (2.4)–(2.7) with boundary conditions (2.2), i.e., to solve
Ldu := Adux+ Bduy = Φ in Ω, (2.18)
Rdu := Cdu = 0 on Γ1, (2.19)
where Φ= (ϕ1, ϕ2, ϕ3, −ϕ1− p)t, the vector Φ = (ϕ1, ϕ2, ϕ3, p)tsolves the stress-pressure
problem (2.16)–(2.17), and Ad= 1 0 0 0 0 1 0 0 , Bd= 0 0 1 0 0 0 0 1 , Cd= 1 0 0 1 .
For constructing numerical solutions to these first-order problems, we will apply the least-squares principles in connection with finite element techniques to both stages.
It is interesting to point out that although the displacement system involves four first-order equations (2.4)–(2.7) with two unknown functions u1 and u2, those first-order equations are pairwise dependent, according to (2.10) and (2.11).
We shall require some function spaces defined on Ω throughout this article [28, 29]. We let Hs(Ω), s ≥ 0 integer, denote the Sobolev space of functions that have square-integrable
derivatives of order up to s on Ω; as usual, L2(Ω) := H0(Ω). The associated inner product and norm are given by
(u, v)s= X |α|≤s Z Ω∂ αu · ∂αv, kuks= p (u, u)s,
respectively. For the product space [Hs(Ω)]m, the corresponding usual inner product and norm
are also denoted by (·, ·)sand k · ks, respectively, when there is no chance for confusion. Let Hs
0(Ω) be the closure of D(Ω) in Hs(Ω), where D(Ω) denotes the linear space of infinitely differentiable functions on Ω with compact support. We denote by H−s(Ω) the dual space of Hs 0(Ω) normed by kuk−s= sup 06=v∈Hs 0(Ω) hu, vi kvks,
where h·, ·i denotes the duality pairing.
The existence, uniqueness, and smoothness of the solution of original second-order problem (2.1)–(2.3) with smooth data are well-known (see, e.g., [30, 31]). Thus, it is reasonable to assume that the stress-pressure problem (2.16)–(2.17) and the displacement problem (2.18)–(2.19) have unique (strong) solutions Φ ∈ [H1(Ω)]4 and u ∈ [H1(Ω)]2, respectively, for given functions
F ∈ [L2(Ω)]4and G ∈ [L2(∂Ω)]2. It is also understood that, when = 0, we further require R
Ωp = 0 here, as well as in the approximations. Actually, the unique solvability of problem (2.18)–(2.19) can be ensured by virtue of Eqs. (2.4)–(2.7) and the boundary requirement u = 0 on Γ1, provided measure (Γ1) > 0.
For simplicity, we shall also assume that the boundary data G in (2.17) is identical to 0, i.e., g = (g1, g2)t= 0 on Γ2. This can be achieved under some suitable assumptions. For example, assume there exist ψ1, ψ2∈ H1(Ω) such that the traces of ψ1and ψ2on Γ are given, respectively, by ψ1= 0 on Γ1; g1n1− g2n2 on Γ2, and ψ2= 0 on Γ1; g2n1+ g1n2 on Γ2. Define ϕ∗ 1= ψ2µ1, ϕ∗2=ψµ2. By the change of variables ˜Φ = ( ˜ϕ1, ˜ϕ2, ˜ϕ3, ˜p)t, where
˜
the stress-pressure problem (2.16)–(2.17) can be transformed into the following form:
Lsp˜Φ := Asp˜Φx+ Bsp˜Φy= ˜F in Ω, Rsp˜Φ := Csp˜Φ = 0 on Γ,
which is the desired result.
III. TWO-STAGE LEAST-SQUARES PROCEDURE We first define two function spaces for our problems,
S = {Ψ ∈ [H1(Ω)]4, R
spΨ = 0 on Γ}, (3.1)
V = {v ∈ [H1(Ω)]2; R
dv = 0 on Γ1}, (3.2)
and then define the least-squares quadratic functional Jsp: S → R by
Jsp(Ψ) = kLspΨ − F k20= k(AspΨx+ BspΨy) − F k20 ∀Ψ ∈ S. (3.3) It is evident that the exact solution Φ ∈ S of the stress-pressure problem (2.16)–(2.17) minimizes (3.3), since Jsp(Φ) = 0, and a zero minimizer of the functional Jsp on S solves problem
(2.16)–(2.17). Thus, the least-squares method for (2.16)–(2.17) is defined to be the following minimization problem:
Seek Φ ∈ S such that Jsp(Φ) = minΨ∈SJsp(Ψ). (3.4)
Taking the first variation, we can find that problem (3.4) is equivalent to
Seek Φ ∈ S such that Bsp(Φ, Ψ) = Fsp(Ψ) ∀Ψ ∈ S, (3.5)
where Bsp(Φ, Ψ) = Z Ω(AspΦx+ BspΦy) · (AspΨx+ BspΨy) ∀Φ, Ψ ∈ S, (3.6) Fsp(Ψ) = Z ΩF · (AspΨx+ BspΨy) ∀Ψ ∈ S. (3.7) Therefore, a least-squares finite element approximation to the solution of problem (2.16)– (2.17) is defined by
Seek Φh∈ Shrsuch that Bsp(Φh, Ψh) = Fsp(Ψh) ∀Ψh∈ Shr, (3.8)
where the finite-dimensional subspace Sr
h⊂ S (with r ≥ 0) is assumed to possess the following
approximation property: for every Ψ ∈ S ∩ [Hr+1(Ω)]4, there exists Ψh∈ Sr
hsuch that kΨ − Ψhk0+ hkΨ − Ψhk1≤ Chr+1kΨkr+1, (3.9)
where C is a positive constant independent of Ψ and h. In what follows, C will denote a positive constant always independent of h, not necessarily the same in different occurrences.
After the stress problem (2.16)–(2.17) is solved by using the least-squares finite element scheme (3.8), our second stage is to solve the displacement problem (2.18)–(2.19) approximately. Define the following least-squares functional
Jd(v) = kLdv − Φk20= k(Advx+ Bdvy) − Φk20 ∀v ∈ V, (3.10) where Φ= (ϕ
1, ϕ2, ϕ3, −ϕ1− p)t, with Φ = (ϕ1, ϕ2, ϕ3, p)tbeing the solution of problem (2.16)–(2.17). Similar to the least-squares method for the stress-pressure problem, we define the following minimization problem:
Seek u ∈ V such that Jd(u) = min
v∈VJd(v), (3.11)
or, equivalently,
Seek u ∈ V such that Bd(u, v) = Fd(v) ∀v ∈ V, (3.12)
where Bd(u, v) = Z Ω(Adux+ Bduy) · (Advx+ Bdvy) ∀u, v ∈ V, (3.13) Fd(v) = Z ΩΦ · (Adv x+ Bdvy) ∀v ∈ V. (3.14)
Since the data function Φ can be obtained only through the numerical scheme (3.8), the
associated least-squares approximate scheme for (2.18)–(2.19) is defined by
Seek uh∈ Vhpsuch that Bd(uh, vh) = ˜Fd(vh) ∀vh∈ Vhp, (3.15)
where ˜ Fd(vh) = Z ΩΦ h· (Advhx+ Bdvhy) ∀vh∈ Vhp, (3.16) Φ
h = (ϕ1h, ϕ2h, ϕ3h, −ϕ1h− ph)t, the vector Φh = (ϕ1h, ϕ2h, ϕ3h, ph)tis the solution of problem (3.8), and the finite-dimensional subspace Vhp ⊂ V (for p ≥ 0) is also assumed to be
equipped with the following approximation property: for every v ∈ V ∩[Hp+1(Ω)]2, there exists vh∈ Vhpsuch that
kv − vhk0+ hkv − vhk1≤ Chp+1kvkp+1, (3.17)
where C is a positive constant independent of v and h.
It is easily seen that Bsp(·, ·) and Bd(·, ·) define two inner products on S × S and V × V,
respectively. The positive-definiteness is ensured by the fact that the stress-pressure problem (2.16)–(2.17) and the displacement problem (2.18)–(2.19) possess a unique solution for each given smooth function F ∈ [L2(Ω)]4. Denote the reduced norms, respectively, by
kΨksp=qBsp(Ψ, Ψ) ∀Ψ ∈ S, (3.18)
kvkd=
p
Bd(v, v) ∀v ∈ V. (3.19)
Then, evidently, there exists a positive constant C such that
kvkd≤ Ckvk1 ∀v ∈ V, (3.21) since both Lspand Ldare first-order differential operators with constant coefficients.
We have the following fundamental properties of the first stage (3.8).
Theorem 3.1. Let Φ ∈ S be the solution of the stress-pressure problem (2.16)–(2.17) with the given functions F ∈ [L2(Ω)]4and G = 0.
(i) Problem (3.8) has a unique solution Φh∈ Shrsatisfying the following stability estimate:
kΦhksp≤ kF k0. (3.22)
(ii) The matrix of the linear system associated with problem (3.8) is symmetric and positive
definite.
(iii) The following orthogonality relation holds:
Bsp(Φ − Φh, Ψh) = 0 ∀Ψh∈ Shr. (3.23)
(iv) The approximate solution Φhis a best approximation of Φ in the k · ksp-norm, that is, kΦ − Φhksp= inf
Ψh∈ShrkΦ − Ψhksp. (3.24)
(v) If Φ ∈ S ∩ [Hr+1(Ω)]4, then
kLspΦh− F k0= kΦ − Φhksp≤ ChrkΦkr+1. (3.25)
Proof. To prove the unique solvability, it suffices to prove the uniqueness of the solution, since
Sr
his a finite-dimensional space. Let Φh be a solution of (3.8). Then by the Cauchy–Schwarz
inequality, we have
kΦhk2sp = Bsp(Φh, Φh) = (F, LspΦh)0
≤ kF k0kLspΦhk0 = kF k0kΦhksp,
which implies (3.22). Consequently, the solution Φhof problem (3.8) is unique.
Assertion (ii) follows from the fact that the inner product Bsp(·, ·) is symmetric and positive
definite. Subtracting the equation in (3.8) from the equation in (3.5), we get (3.23). To prove (iv), by (3.23) and the Cauchy–Schwarz inequality, we obtain
kΦ − Φhk2sp = Bsp(Φ − Φh, Φ − Φh)
= Bsp(Φ − Φh, Φ − Ψh)
≤ kΦ − ΦhkspkΦ − Ψhksp ∀Ψh∈ Sr h,
which implies (3.24).
Finally, let Ψh ∈ Shrsuch that (3.9) holds with Ψ replaced by Φ. Together with (3.24) and
(3.20), we can obtain (3.25).
Similar to Theorem 3.1 with minor modifications, we have the following results for second stage (3.15).
(i) Problem (3.15) has a unique solution uh∈ Vhpsatisfying the following stability estimate: kuhkd≤ kΦ
hk0, (3.26)
where Φ
h= (ϕ1h, ϕ2h, ϕ3h, −ϕ1h−ph)t, and Φh= (ϕ1h, ϕ2h, ϕ3h, ph)tis the solution
of problem (3.8).
(ii) The matrix of the linear system associated with problem (3.15) is symmetric and positive
definite.
(iii) The following relation holds:
Bd(u − uh, vh) = (Φ− Φ
h, Ldvh)0 ∀vh∈ Vhp, (3.27) where Φ= (ϕ
1, ϕ2, ϕ3, −ϕ1− p)t, and Φ = (ϕ1, ϕ2, ϕ3, p)tis the solution of problem (2.16)–(2.17).
In the following two sections, we will prove that the k · ksp-norm and the k · kd-norm are
equivalent to the k · k1-norm in the respective spaces S and V. By (3.22) and (3.26), we will have the following corollary.
Corollary 3.1. The two-stage least-squares finite element scheme (3.8) and (3.15) is stable with respect to the k · k1-norm, i.e.,
kΦhk1≤ CkF k0, (3.28)
kuhk1≤ CkF k0, (3.29)
where C is a positive constant independent of h.
IV. A PRIORI ESTIMATES
In this section we shall apply the theory given by Wendland [26] to derive coercive type a
priori estimates for the solution Φ to the stress-pressure problem (2.16)–(2.17). Following these
estimates, the error estimates for the least-squares finite element approximation (3.8) can be obtained.
We shall show that Lspis an elliptic operator in the sense of Petrovski, and that the boundary
operator Rspin (2.17) satisfies the Lopatinski condition. Then (Lsp, Rsp) is a regular elliptic
system and so, by [26], it is a Fredholm operator with zero nullity, which enables us to get the a
priori estimates (cf. Theorem 4.1 below).
For all (ξ, η) ∈ R2and (ξ, η) 6= (0, 0), we have
det(ξAsp+ ηBsp) = 2µ2(1 + )(ξ2+ η2)26= 0.
Thus, (2.16) is an elliptic system in the sense of Petrovski. Taking (ξ, η) = (1, 0), we find that
Aspis nonsingular and
A−1 sp = − 2µ(1+) 0 0 1+1 0 0 −1 0 0 −1 µ 1 0 1 2µ(1+) 0 0 1+1 .
Then the original elliptic system (2.16) can be transformed into the following form: Φx+ ˜BspΦy = ˜F in Ω, where ˜ Bsp= A−1spBsp= 0 2(1+) 2(1+)2+ 0 −1 0 0 0 −1 0 0 −2(1 + ) 0 − 1 2(1+) 2(1+)1 0 , ˜ F = A−1 spF = − 2µ(1+)f1 0 −1 µf2 1 2µ(1+)f1 .
We now check the Lopatinski condition as follows: after elementary operations, we can find that the eigenvalues of matrix ˜Bt
spare i and −i, each with multiplicity two. Consider the eigenvalue τ+ = i in the upper half-plane, to which there is a pair of linearly independent generalized eigenvectors p1and p2of ˜Bspt obtained from
˜ Bt spp1− τ+p1= 0, ˜ Bt spp2− τ+p2= p1, where p1= (0, 1, −1, −2(1 + )i)t, p2= −4(1 + )2 + , 0,4(1 + )2 + i, −2(1 + )(2 + 3)2 + t . Notice that P = (p1, ¯p1, p2, ¯p2)t is nonsingular. Let Q = (q1, ¯q1, q2, ¯q2) be the inverse matrix of P. Then
Q = − 2+3 8(1+)i 8(1+)2+3 i −8(1+)2+ −8(1+)2+ 1 2 12 −8(1+)2+ i 8(1+)2+ i 0 0 − 2+ 8(1+)i 8(1+)2+ i 1 4(1+)i −4(1+)1 i 0 0 .
Now we have det{2Csp(q1, q2)} = (2+)(2+3) 16(1+)2 (n1+ n2i)26= 0 on Γ1; −16(1+)(2+)22(n1+ n2i)26= 0 on Γ2,
for all ≥ 0 and (n1, n2) 6= (0, 0). That is, the Lopatinski condition is fulfilled for our boundary conditions (2.17). Thus, we have the following theorem.
Theorem 4.1. For the boundary value problem (2.16)–(2.17), we have the following a priori estimates: for each l ≥ 0 there is a constant C > 0 such that if Ψ ∈ [Hl+1(Ω)]4, then
kΨkl+1≤ C{kLspΨkl+ kRspΨkl+1
2}. (4.1)
By an interpolation argument in [35] (see also [26, Lemma 8.2.1]), the inequalities (4.1) can be extended to the case l ≥ −1. Taking l = 1, l = 0, and l = −1 in (4.1), we have, respectively,
kΨk2≤ CkLspΨk1 ∀Ψ ∈ S ∩ [H2(Ω)]4, (4.2)
kΨk1≤ CkLspΨk0 ∀Ψ ∈ S, (4.3)
kΨk0≤ CkLspΨk−1 ∀Ψ ∈ S. (4.4)
The a priori estimates (4.2), (4.3), and (4.4) play crucial roles for the least-squares error estimates of the stresses and pressure in the next section.
Remark 4.1. It is unclear whether the constant C in (4.2)–(4.4) is independent of the non-negative parameter , since the constant in (4.1) is not explicitly known (cf. [32], page 74, Remark 2).
V. ERROR ESTIMATES
It is easily seen that the bilinear form Bspis continuous on S × S; the coercivity of Bspfollows
from (4.3). Thus, we first obtain the following result.
Theorem 5.1. Let Φ ∈ S and Φh ∈ Shrbe the solutions of problems (2.16)–(2.17) and (3.8), respectively. Assume that Φ ∈ [Hr+1(Ω)]4. Then there exists a positive constant C independent
of h and Φ such that
kΦ − Φhk1≤ ChrkΦkr+1. (5.1)
Proof. Since the k · ksp-norm is equivalent to the k · k1-norm on the space S, the assertion follows from (3.25) immediately.
For deriving the optimal L2-estimates, we need the following regularity assumption that we shall use in the subsequent result.
Assumption (A1). For any Ψ ∈ [H1
0(Ω)]4, the unique solution ˜Φ of the problem
Lsp˜Φ = Ψ in Ω,
belongs to S ∩ [H2(Ω)]4.
Evidently, this is a reasonable assumption because the differential operator Lspis of first order
and the data function Ψ is in [H1 0(Ω)]4. Theorem 5.2. Let Φ ∈ S and Φh ∈ Sr
hbe the solutions of problems (2.16)–(2.17) and (3.8), respectively. Assume that Φ ∈ [Hr+1(Ω)]4and that assumption (A1) holds. Then there exists a
positive constant C independent of h and Φ such that
kΦ − Φhk0≤ Chr+1kΦkr+1. (5.3)
Proof. Let Ψ ∈ [H1
0(Ω)]4and let ˜Φ ∈ S ∩[H2(Ω)]4be the corresponding solution to problem (5.2). Then |(Lsp(Φ − Φh), Ψ)0| = |(Lsp(Φ − Φh), Lsp˜Φ)0| = |(Lsp(Φ − Φh), Lsp(˜Φ − Ψh))0| ∀Ψh∈ Shr (by (3.23)) ≤ kLsp(Φ − Φh)k0kLsp(˜Φ − Ψh)k0 ∀Ψh∈ Shr ≤ CkΦ − Φhk1k˜Φ − Ψhk1 ∀Ψh∈ Shr ≤ ChkΦ − Φhk1k˜Φk2 (by (3.9)) ≤ ChkΦ − Φhk1kLsp˜Φk1 (by (4.2)) = ChkΦ − Φhk1kΨk1.
In addition, the L2 inner product (L
sp(Φ − Φh), Ψ)0 defines a bounded linear functional on [H1
0(Ω)]4, since
|(Lsp(Φ − Φh), Ψ)0| ≤ kLsp(Φ − Φh)k0kΨk1 ∀Ψ ∈ [H01(Ω)]4. Therefore,
kLsp(Φ − Φh)k−1≤ ChkΦ − Φhk1. (5.4) Combining (5.4) with (4.4) and (5.1), we can readily conclude estimate (5.3).
The results of Theorem 5.1 and Theorem 5.2 indicate that the rates of convergence for the stresses and pressure are optimal, both in the H1-norm and in the L2-norm. We now estimate the rates of convergence of the approximations for the displacements. The continuity of the bilinear form Bdcan be obtained easily. For any v = (v1, v2)t ∈ V = {v ∈ [H1(Ω)]2; v = 0 on Γ1}, we have Bd(v, v) = Z Ω ∂v1 ∂x 2 + ∂v1 ∂y 2 + ∂v2 ∂x 2 + ∂v2 ∂y 2 = k∇vk2 0.
It follows from the Poincar´e–Friedrichs inequality that
Bd(v, v) ≥ Ckvk2
1, (5.5)
i.e., Bdis coercive on V × V. Similar to Theorem 5.1, we have the optimal order of convergence
in the H1-norm for the displacements.
Theorem 5.3. Let Φ ∈ S, u ∈ V, and uh ∈ Vhpbe the solutions of problems (2.16)–(2.17),
(2.18)–(2.19), and (3.15), respectively. Assume Φ ∈ [Hr+1(Ω)]4 and u ∈ [Hp+1(Ω)]2. Then
there exists a positive constant C independent of h, u, and Φ such that
Proof. For any vh∈ Vhp, by (5.5) and (3.27) we have Ckuh− vhk21 ≤ Bd(uh− vh, uh− vh) = Bd(u − vh, uh− vh) − Bd(u − uh, uh− vh) = Bd(u − vh, uh− vh) − (Φ− Φh, Ld(uh− vh))0 ≤ Cku − vhk1kuh− vhk1+ CkΦ− Φhk0kuh− vhk1 ≤ C(ku − vhk1+ kΦ − Φhk0)kuh− vhk1, which implies kuh− vhk1≤ C(ku − vhk1+ kΦ − Φhk0) ∀vh∈ Vhp. Thus, ku − uhk1 ≤ ku − vhk1+ kuh− vhk1 ≤ C(ku − vhk1+ kΦ − Φhk0) ∀vh∈ Vhp.
Choose vh∈ Vhpsuch that (3.17) holds with v replaced by u. Then we obtain ku − uhk1 ≤ C(ku − vhk1+ kΦ − Φhk1)
≤ C(hpkukp+1+ hrkΦkr+1).
This completes the proof.
Similar to the derivation of Theorem 5.2, we shall use the Aubin–Nitsche trick [28, 33] to establish the optimal L2-estimates for the displacements. For each h, consider the following adjoint problem:
Find ˜u ∈ V such that Bd(˜u, v) = (u − uh, v)0 ∀v ∈ V. (5.7) Note that the right-hand side of the equation in (5.7) defines a bounded linear functional on V. Thus, the unique solvability of problem (5.7) is ensured by the Lax–Milgram lemma, since Bdis
coercive on V × V. We now assume the following regularity assumption (cf. [14]).
Assumption (A2). Assume the unique solution ˜u of problem (5.7) belongs to [H2(Ω)]2∩ V
and there exists a positive constant C independent of ˜u and u − uhsuch that
k˜uk2≤ Cku − uhk0. (5.8)
Then we have the following optimal L2-estimates for the displacements.
Theorem 5.4. Let Φ ∈ S, u ∈ V, and uh ∈ Vhpbe the solutions of problems (2.16)–(2.17),
(2.18)–(2.19), and (3.15), respectively. Assume that Φ ∈ [Hr+1(Ω)]4, u ∈ [Hp+1(Ω)]2, and
regularity assumptions (A1) and (A2) hold. Then there exists a positive constant C independent of h, u, and Φ such that
ku − uhk0≤ C(hp+1kukp+1+ hr+1kΦkr+1). (5.9) Proof. Choosing v = u − uh∈ V in (5.7), we have
Bd(˜u, u − uh) = ku − uhk20, which together with (3.27) enables us to obtain
Bd(˜u − wh, u − uh) = ku − uhk2
It follows that
ku − uhk20 = Bd(˜u − wh, u − uh) + (Φ− Φh, Ldwh)0
= Bd(˜u − wh, u − uh) + (Φ− Φh, Ld(wh− ˜u))0+ (Φ− Φh, Ld˜u)0 ≤ k˜u − whk1ku − uhk1+ kΦ− Φhk0k˜u − whk1+ kΦ− Φhk0k˜uk2, for all wh∈ Vhp. Choose wh∈ Vhpso that
k˜u − whk1≤ Chk˜uk2. Hence, together with (5.8) we have
ku − uhk2
0≤ Chku − uhk0ku − uhk1+ ChkΦ − Φhk0ku − uhk0+ CkΦ − Φhk0ku − uhk0, which implies
ku − uhk0≤ C(hku − uhk1+ hkΦ − Φhk0+ kΦ − Φhk0). This completes the proof.
VI. CONDITION NUMBERS
In this section, we shall give estimates for the condition numbers of the linear systems arising from problem (3.8) and problem (3.15). Recall that the condition number for a symmetric and positive definite m × m matrix M is defined by
condition number of M = λmax λmin =
max ρ(Ξ) min ρ (Ξ),
where λmaxand λminare the largest and smallest eigenvalues of M, and ρ(Ξ) is the Rayleigh quotient,
ρ(Ξ) := ΞΞtMΞtΞ ∀Ξ = (ξ1, . . . , ξm)t∈ Rm, Ξ 6= 0.
We shall assume that the respective bases {Φ1, . . . , ΦK} and {u1, . . . , uk} of the finite element spaces Sr
hand Vhpare chosen so that the following conditions hold: there exist positive constants
Λ1and Λ2such that for all ξ1, . . . , ξK, η1, . . . , ηk ∈ R, Λ1h2 K X i=1 ξ2 i ≤ K X i=1 ξiΦi 2 0 ≤ Λ2h2 K X i=1 ξ2 i, (6.1) Λ1h2 k X i=1 η2 i ≤ k X i=1 ηiui 2 0 ≤ Λ1h2 k X i=1 η2 i. (6.2)
The above conditions hold for most finite element spaces Sr
hand Vhp. If, in addition, the
corre-sponding regular family {Th} of triangulations of ¯Ω is quasi-uniform [28, 34], i.e., there exists a
positive constant C independent of h such that
h ≤ C diam(Ωh
TABLE I. The approximations Φhwith E = 2.5 and ν = 0.25 ( = 2.0).
1/h L2-error RelErr L2-rate k · k
sp-error RelErr k · ksp-rate
2 0.85600 2.91285e-1 — 7.15156 4.58279e-1 —
4 0.24756 8.42423e-2 1.79 3.52245 2.25722e-1 1.02
8 0.06657 2.26540e-2 1.89 1.76541 1.13130e-1 1.00
16 0.01708 5.81339e-3 1.96 0.88398 5.66466e-2 1.00
32 0.00431 1.46587e-3 1.99 0.44219 2.83361e-2 1.00
then we have the following inverse estimates: K X i=1 ξiΦi 2 1 ≤ Ch−2 K X i=1 ξiΦi 2 0 ≤ CΛ2 K X i=1 ξ2 i, (6.3) k X i=1 ηiui 2 1 ≤ Ch−2 k X i=1 ηiui 2 0 ≤ CΛ2 k X i=1 η2 i, (6.4)
where C is a positive constant independent of h.
Theorem 6.1. Under conditions (6.1) and (6.3) [respectively, (6.2) and (6.4)] the condition number of the linear system arising from problem (3.8) [respectively, problem (3.15)] is O(h−2).
Proof. Let Φh:=PKi=1ξiΦi ∈ Shr. Since the bilinear form Bsp(·, ·) is coercive on S × S,
by (6.1) we have Bsp(Φh, Φh) ≥ CkΦhk2 1≥ CkΦhk20≥ CΛ1h2 K X i=1 ξ2 i.
On the other hand, by the continuity of Bsp(·, ·) on S × S, we get from (6.3) that Bsp(Φh, Φh) ≤ CkΦhk21≤ CΛ2 K X i=1 ξ2 i.
Thus, λmax ≤ CΛ2 and λmin ≥ CΛ1h2, and so the condition number for problem (3.8) is O(h−2). The estimates of the condition number for problem (3.15) can be achieved in a
similar way.
TABLE II. The approximations uhwith E = 2.5 and ν = 0.25 ( = 2.0).
1/h L2-error RelErr L2-rate k · k
d-error RelErr k · kd-rate
2 0.17245 2.43881e-1 — 1.40921 4.48564e-1 —
4 0.04299 6.08011e-2 2.00 0.70905 2.25697e-1 0.99
8 0.01075 1.52010e-2 2.00 0.35569 1.13221e-1 1.00
16 0.00269 3.79905e-3 2.00 0.17801 5.66632e-2 1.00
TABLE III. Rates of convergence in the k · ksp-norm with E = 2.5 and small .
ν = 0.49 ν = 0.499 ν = 0.4999 ν = 0.49999 ν = 0.499999
1/h ' 4.1e-2 ' 4.0e-3 ' 4.0e-4 ' 4.0e-5 ' 4.0e-6
2 — — — — —
4 0.97 0.96 0.96 0.96 0.96
8 0.98 0.98 0.97 0.97 0.97
16 0.99 0.99 0.99 0.99 0.99
32 1.00 1.00 1.00 1.00 1.00
VII. NUMERICAL EXPERIMENTS
We shall present a simple example solved by using our two-stage least-squares finite element scheme (3.8) and (3.15). To simplify the numerical implementation, we shall assume that Ω = (0, 1) × (0, 1), Γ1 = ∂Ω, and the square domain Ω is uniformly partitioned into a set of 1/h2 square subdomains Ωh
i with side-length h. The problem we present has the smooth exact solution,
ϕ1 ϕ2 ϕ3 p u1 u2 = π cos(πx) sin(πy) π sin(πx) cos(πy) π cos(πx) sin(πy) −π
(cos(πx) sin(πy) + sin(πx) cos(πy))
sin(πx) sin(πy) sin(πx) sin(πy) . (7.1)
Substituting (7.1) into (2.16)–(2.17), we have G = (0, 0)t, F = (f
1, f2, 0, 0)t, and f1= 2µπ2 3 2 + 1 sin(πx) sin(πy) − 1 2 + 1 cos(πx) cos(πy) , f2= 2µπ2 3 2 + 1 sin(πx) sin(πy) − 1 2 + 1 cos(πx) cos(πy) .
Piecewise bilinear finite elements are applied for all the unknowns. For the case of Poisson’s ratio ν = 0.25 and Young’s modulus E = 2.5, the results are collected in Table I and Table II, where RelErr denotes the relative error and, for simplicity, the data function Φ
his replaced by the
exact function Φin the second stage (3.15). Since the k · k
sp-norm and the k · kd-norm are both
equivalent to the H1-norm on the spaces S and V, respectively, the numerical results in Table I and Table II indicate that the two-stage least-squares procedure (3.8) (3.15) achieves optimal convergence both in the L2-norm and in the H1-norm for all the unknowns.
TABLE IV. Rates of convergence in the L2-norm with E = 2.5 and small .
ν = 0.49 ν = 0.499 ν = 0.4999 ν = 0.49999 ν = 0.499999
1/h ' 4.1e-2 ' 4.0e-3 ' 4.0e-4 ' 4.0e-5 ' 4.0e-6
2 — — — — —
4 1.78 1.76 1.76 1.76 1.76
8 1.89 1.89 1.89 1.89 1.89
16 1.96 1.95 1.95 1.95 1.95
The behavior of convergence for the stresses and pressure influenced by the nonnegative parameter is particularly examined. Table III and Table IV exhibit that, except on very coarse meshes, the optimal convergence is still essentially ensured for various values of the parameter, even for nearly incompressible elasticity. That is, computational results in Table III and Table IV indicate that the behavior of convergence is uniform in the nonnegative parameter.
The authors would like to thank two anonymous referees for their careful reading of the manuscript and their fruitful comments and suggestions.
References
1. A. K. Aziz, R. B. Kellogg, and A. B. Stephens,‘‘Least squares methods for elliptic systems,’’Math. Comp. 44, 53–70 (1985).
2. P. B. Bochev and M. D. Gunzburger,‘‘Analysis of least squares finite element methods for the Stokes equations,’’Math. Comp. 63, 479–506 (1994).
3. P. B. Bochev and M. D. Gunzburger,‘‘Least-squares methods for the velocity-pressure-stress formu-lation of the Stokes equations,’’Comput. Methods Appl. Mech. Engrg. 126, 267–287 (1995).
4. J. H. Bramble and J. E. Pasciak,‘‘Least-squares methods for Stokes equations based on a discrete minus one inner product,’’J. Comp. Appl. Math. 74, 155–173 (1996).
5. C. L. Chang,‘‘A mixed finite element method for the Stokes problem: an acceleration-pressure for-mulation,’’Appl. Math. Comp. 36, 135–146 (1990).
6. C. L. Chang and J. J. Nelson,‘‘Least-squares finite element method for the Stokes problem with zero residual of mass conservation,’’SIAM J. Numer. Anal. 34, 480–489 (1997).
7. C. L. Chang and S.-Y. Yang,‘‘Piecewise linear approximations to the Stokes problem with velocity boundary conditions,’’preprint, November 1996, to appear.
8. C. L. Chang, S.-Y. Yang, and C.-H. Hsu,‘‘A least-squares finite element method for incompressible flow in stress-velocity-pressure version,’’Comput. Methods Appl. Mech. Engrg. 128, 1–9 (1995).
9. T.-F. Chen,‘‘On least squares approximations to compressible flow problems,’’Numerical Methods for PDEs 12, 207–228 (1986).
10. T.-F. Chen and G. J. Fix,‘‘Least squares finite element simulation of transonic flows,’’Appl. Numer. Math. 2, 399–408 (1986).
11. B.-N. Jiang and C. L. Chang,‘‘Least-squares finite elements for the Stokes problem,’’Comput. Methods Appl. Mech. Engrg. 78, 297–311 (1990).
12. B.-N. Jiang, T. L. Lin, and L. A. Povinelli,‘‘Large-scale computation of incompressible viscous flow by least-squares finite element method,’’Comput. Methods Appl. Mech. Engrg. 114, 213–231 (1994).
13. C. L. Chang and M. D. Gunzburger, ‘‘A subdomain-Galerkin/least squares method for first-order elliptic systems in the plane,’’SIAM J. Numer. Anal. 27, 1197–1211 (1990).
14. G. J. Fix and M. E. Rose,‘‘A comparative study of finite element and finite difference methods for Cauchy–Riemann type equations,’’SIAM J. Numer. Anal. 22, 250–261 (1985).
15. L. P. Franca and T. J. R. Hughes,‘‘Two classes of mixed finite element methods,’’Comput. Methods Appl. Mech. Engrg. 69, 89–129 (1988).
16. L. P. Franca and R. Stenberg,‘‘Error analysis of some Galerkin least squares methods for the elasticity equations,’’SIAM J. Numer. Anal. 28, 1680–1697 (1991).
17. S.-Y. Yang and J.-L. Liu,‘‘Analysis of least squares finite element methods for a parameter-dependent first-order system,’’preprint, March 1997, to appear.
18. C. L. Chang,‘‘Finite element method for the solution of Maxwell’s equations in multiple media,’’ Appl. Math. Comp. 25, 89–99 (1988).
19. J. S. Li, Z. Y. Yu, X. Q. Xiang, W. P. Ni, and C. L. Chang,‘‘Least-squares finite element method for electromagnetic fields in 2-D,’’Appl. Math. Comp. 58, 143–167 (1993).
20. D. M. Bedivan and G. J. Fix,‘‘Least squares methods for optimal shape design problems,’’Computer Math. Applic. 30, 17–25 (1995).
21. J. H. Bramble, R. D. Lazarov, and J. E. Pasciak,‘‘A least-squares approach based on a discrete minus one inner product for first order systems,’’Math. Comp. 66, 935–955 (1997).
22. Z. Cai, R. Lazarov, T. A. Manteuffel, and S. F. McCormick,‘‘First-order system least squares for second-order partial differential equations: Part I,’’SIAM J. Numer. Anal. 31, 1785–1799 (1994).
23. A. I. Pehlivanov, G. F. Carey, and R. D. Lazarov,‘‘Least-squares mixed finite elements for second-order elliptic problems,’’SIAM J. Numer. Anal. 31, 1368–1377 (1994).
24. F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer–Verlag, New York, 1991. 25. J. Pitk¨aranta and R. Stenberg,‘‘Analysis of some mixed finite element methods for plane elasticity
equations,’’Math. Comp. 41, 399–423 (1983).
26. W. L. Wendland, Elliptic Systems in the Plane, Pitman, London, 1979.
27. A. K. Aziz and A. Werschulz,‘‘On the numerical solution of Helmholtz’s equation by the finite element method,’’SIAM J. Numer. Anal. 17, 681–686 (1980).
28. P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978. 29. J. L. Lions and E. Magenes, Nonhomogeneous Elliptic Boundary Value Problems and Applications,
Vol. I, Springer–Verlag, Berlin, 1972.
30. G. Fichera, Linear Elliptic Differential Systems and Eigenvalue Problems, Lecture Notes in Mathematics 8, Springer–Verlag, Berlin, 1965.
31. G. Fichera,‘‘Existence theorems in elasticity,’’Handbuch der Physik 6, a/2, 347–389 (1972).
32. S. Agmon, A. Douglis, and L. Nirenberg,‘‘Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions II,’’Comm. Pure Appl. Math. 17, 35–92
(1964).
33. G. Strang and G. J. Fix, An Analysis of the Finite Element Method, Prentice–Hall, Englewood Cliffs, NJ, 1973.
34. C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cam-bridge University Press, 1990.
35. A. S. Dikanskij,‘‘Conjugate problems of elliptic differential and pseudo differential boundary value problems in a bounded domain,’’Math. USSR Sbornik 20, 67–83 (1973).
