- Email: [email protected]
On efficient least-squares finite element methods for convection-dominated problems
Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
Contents lists available at ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
On efficient least-squares finite element methods for convection-dominated problems Po-Wen Hsieh, Suh-Yuh Yang * Department of Mathematics, National Central University, Jhongli City 32001, Taiwan
a r t i c l e
i n f o
Article history: Received 2 September 2008 Accepted 30 September 2009 Available online 9 October 2009 MSC: 65N12 65N15 65N30 Keywords: Convection-dominated problems Boundary and interior layers Finite element methods Least-squares Stabilized methods Streamline diffusion
a b s t r a c t This paper focuses on the least-squares finite element method and its three variants for obtaining efficient numerical solutions to convection-dominated convection–diffusion problems. The coercivity estimates for the corresponding homogeneous least-squares energy functionals are derived, and based on which error estimates are established. One of the common advantages of these least-squares methods is that the resulting linear system is symmetric and positive definite. Numerical experiments that demonstrate the theoretical analysis of the developed methods are presented. It was observed that the primitive least-squares method performs poorly for convection-dominated problems while the stabilized, streamline diffusion and negatively stabilized streamline diffusion least-squares methods perform considerably better for interior layer problems, and the negatively stabilized streamline diffusion leastsquares method is able to better capture the boundary layer behavior when compared with other least-squares methods. But all the least-squares methods do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction Let X be a bounded convex polygonal domain in R2 with boundary @ X. We consider the following convection–diffusion problem:
jDu þ a ru ¼ f u¼g
in X;
on @ X;
ð1Þ
where 0 < j 6 1 is the constant viscosity coefficient (or diffusivity), a ¼ ða1 ; a2 Þ> is the velocity field, f is a source function and g is a given boundary data. It is well-known that the solution u of problem (1) can exhibit localized phenomena such as boundary and interior layers, i.e., narrow regions where the derivative of u is very large. This is the case for instance if the problem is convection-dominated, i.e., 0 < j jaj. For simplifying the mathematical analysis, in what follows we assume that a is a constant vector with jaj ¼ 1 and g ¼ 0 on the boundary @ X. It is already known that conventional numerical methods are lacking in either stability or accuracy for the convection-dominated problems. For example, the standard Galerkin method using piecewise linear elements performs poorly for the convection-dominated
* Corresponding author. Tel.: +886 3 4227151x65130; fax: +886 3 4257379. E-mail address: [email protected] (S.-Y. Yang). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.09.029
model, in which large spurious oscillations are frequently exhibited. In order to defeat this difficulty, the subject of stabilized finite element methods has existed for more than twenty years and is still attractive today (see, e.g., [1–3] or recent review [4] by Franca et al.). The stabilized methods are formed by adding to the standard Galerkin method variational terms that are consistent and meshdependent. This additional term not only improves the numerical stability of the Galerkin method but also preserves good accuracy. Another possible alternative is still to apply the standard finite element method to solve the convection-dominated problem by enriching piecewise linear polynomials with bubble functions (see the early tutorial paper [5] by Brezzi and Russo, or [6–8] and many references therein). The bubble function is chosen so that the computed solution satisfies the original second-order differential equation in the interior of each element and vanishes on the boundary of each element. This produces a stabilized-like formulation. However, the resulting linear systems of above-mentioned two approaches are not symmetric. On the other hand, for the past two decades, least-squares finite element methods (LSFEMs) have become more and more frequently used to approximate the solution of first-order system of partial differential equations arising from fluid and solid mechanics. We refer the reader to [9,10] for some of the vast literature. In a broad sense, the least-squares approach is to minimize the
184
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
residual (measured in some norms) of a first-order system of partial differential equations over a suitable function space. In general, such a first-order system is obtained from an originally higher-order system by introducing one or more additional physically unknown variables such as vorticity, stress, or flux (see, e.g., [11–18]). By applying the least-squares principle on this extended first-order system, the smoothness requirement on the function space can then be relaxed. It is already known that the specific features of the least-squares finite element approach which give it advantages relative to, for instance, the mixed finite element approach [19,20], are as follows: it leads to a minimization problem; it is not subject to the Ladyzhenskaya–Babuska–Brezzi condition; simple equal low-order finite elements, such as the continuous linear elements, can be used for the approximation of all unknowns; the resulting linear system is symmetric and positive definite; the value of the homogeneous least-squares functional of the approximate solution provides a practical and sharp a posteriori error estimator at no additional cost, etc. Unfortunately, as we will see in this paper, the primitive LSFEM is not able to achieve good performance for convection-dominated problems. In [21], the authors proposed an exponentially weighted LSFEM to achieve H1 -like accuracy of the discretization and optimal multigrid convergence uniformly in the size of convection. However, their norms incorporate a scaling that has the effect of dramatically changing the influence of the boundary layer. The aim of this paper is to combine the advantages of abovementioned stabilized and least-squares finite element approaches to tackle the convection-dominated problems. We will introduce and analyze the primitive LSFEM and its three variants. To this aim, we have to recast problem (1) into the following first-order system problem by introducing the new dependent variable p ¼ jru in X:
product ðp; qÞdiv :¼ ðp; qÞ0 þ ðr p; r qÞ0 and norm kqk2div :¼ kqk20 þ kr qk20 . In the analysis of the developed methods, we will frequently use the following Green formula:
8 > < r p þ a ru ¼ f p þ jru ¼ 0 > : u¼0
Fððu; pÞ; f Þ ¼ 0 ¼ minfFððv ; qÞ; f Þ : ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞg:
in X; in X;
ð2Þ
on @ X:
We first apply the primitive LSFEM to approximate the solution to problem (2), in which the corresponding least-squares energy functional is defined to be the sum of the squared L2 norms of the residuals of the partial differential equations over an appropriate product space, and the least-squares finite element solution is defined to be the minimizer of the functional over a continuous piecewise polynomial finite element space. As we will see further below, the primitive LSFEM performs poorly for convection-dominated problems. We then introduce three variants of the primitive LSFEM that will be introduced in next several sections, including the stabilized LSFEM, the streamline diffusion LSFEM [22,23] and the negatively stabilized streamline diffusion LSFEM. Numerical simulations will be performed to demonstrate the theoretical analysis of the proposed methods. More importantly, it was observed that the stabilized, streamline diffusion and negatively stabilized streamline diffusion LSFEMs perform considerably better than the primitive LSFEM for interior layer problems, and the negatively stabilized streamline diffusion LSFEM is able to better capture the boundary layer behavior when compared with other least-squares methods. But regrettably, all the least-squares methods do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. Throughout this paper, we will use standard notation and definition for the Sobolev space Hm ðXÞ for non-negative integer m (cf. [19,24]). The standard associated inner product and norm are denoted by ð; Þm and k km , respectively. As usual, L2 ðXÞ ¼ H0 ðXÞ and H10 ðXÞ ¼ fv 2 H1 ðXÞ and v j@ X ¼ 0g. We also need the space Hðdiv; XÞ ¼ fq 2 ½L2 ðXÞ2 and r q 2 L2 ðXÞg with standard inner
ðq; rv Þ0 þ ðr q; v Þ0 ¼ ðq n; v Þ0;@ X
8 q 2 Hðdiv; XÞ and
v 2 H1 ðXÞ; ð3Þ
where n denotes the unit outward normal vector to @ X, and we also need the following Poincaré–Friedrichs inequality [24]: there exists a constant C pf > 0 such that
kv k0 6 C pf krv k0
8
v 2 H10 ðXÞ:
ð4Þ
The remainder of this paper is organized as follows. We introduce and analyze the primitive LSFEM, the stabilized LSFEM, the streamline diffusion LSFEM, and the negatively stabilized streamline diffusion LSFEM in Sections 2–5, respectively. The coercivity estimates for the corresponding homogeneous least-squares functionals are established, and based on which error estimates are derived. Numerical results for various test cases are reported in Section 6. Some concluding remarks are given in Section 7. 2. The primitive LSFEM Define the L2 least-squares energy functional F : H10 ðXÞ Hðdiv; XÞ ! R for the extended first-order problem (2) by
Fððv ; qÞ; f Þ ¼ kr q þ a rv f k20 þ kq þ jrv 0k20 : Note that the energy functional Fð; f Þ is defined to be the sum of the squared L2 -norms of the residuals on the product function space H10 ðXÞ Hðdiv; XÞ. Thus, if ðu; pÞ 2 H10 ðXÞ Hðdiv; XÞ is an exact solution of problem (2), then ðu; pÞ must be a zero minimizer of the functional Fð; f Þ on H10 ðXÞ Hðdiv; XÞ, namely,
Since Fððu; pÞ þ eðv ; qÞ; f Þ is a nonnegative quadratic functional in e 2 R for any given ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ, we have
d Fððu; pÞ þ eðv ; qÞ; f Þje¼0 ¼ 0; de which is equivalent to
Bððu; pÞ; ðv ; qÞÞ ¼ Lððv ; qÞÞ 8 ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ:
ð5Þ
Here the bilinear form Bð; Þ and the linear form LðÞ are defined as following, respectively:
Bððu; pÞ; ðv ; qÞÞ ¼ ðr p þ a ru; r q þ a rv Þ0 þ ðp þ jru; q þ jrv Þ0 ; Lððv ; qÞÞ ¼ ðf ; r q þ a rv Þ0 : Now we consider the finite element formulation of the leastsquares method (5) for problem (2). We introduce the finite element spaces Vh # H10 ðXÞ and Qh # Hðdiv; XÞ which consist of piecewise polynomials of degree less or equal than r (r P 1 integer) over a regular finite element partition Th of the domain X. We denote by hT the diameter of element T 2 Th and set h :¼ maxT2Th hT . The primitive LSFEM for problem (2) is then defined to be the following:
Determine ðuh ; ph Þ 2 Vh Qh such that Bððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Lððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Qh :
ð6Þ
Once the basis functions of the finite element space Vh Qh are chosen, problem (6) is equivalent to solve a linear system problem, An ¼ b, where unknown vector n consists the coordinates of the finite element solution ðuh ; ph Þ with respect to the chosen basis functions. The unique solvability of problem (6) and error estimates of the approximate solution ðuh ; ph Þ is mainly depending on the coercivity
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
(i.e., stability) and boundedness estimates of bilinear form B. Note that Bððv ; qÞ; ðv ; qÞÞ ¼ Fððv ; qÞ; 0Þ for all ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ. Theorem 1. There exist constants C 1 ; C 2 > 0 independent of such that
C1
j2 kv k21 þ j2 kqk2div
j and h
6 Fððv ; qÞ; 0Þ 6 C 2 kv k21 þ kqk2div ;
According to the orthogonality property (13), we have
Bððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ ¼ Bððuh u; ph pÞ þ ðu uI ; p pI Þ; ðuh uI ; ph pI ÞÞ ¼ Bððu uI ; p pI Þ; ðuh uI ; ph pI ÞÞ
ð7Þ
1
6 B2 ððu uI ; p pI Þ; ðu uI ; p pI ÞÞ 1
for all ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ.
B2 ððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ;
Proof. The boundedness estimate can be easily proven. We proceed to show the validity of the coercivity estimate, namely, the left hand side of (7). Let a be a positive constant that will be determined later. Applying the Poincaré–Friedrichs inequality (4), we have
Fððv ; qÞ; 0Þ ¼ kr q þ a rv av k20 þ 2aðr q þ a rv ; v Þ0 a2 kv k20 þ kq þ jrv arv k20 þ 2aðq þ jrv ; rv Þ0
1
B2 ððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ 1
6 B2 ððu uI ; p pI Þ; ðu uI ; p pI ÞÞ: From above inequality, Theorem 1, interpolation property (12) and the boundedness of the bilinear form B on the space H10 ðXÞ Hðdiv; XÞ,
jkuh uI k1 þ jkph pI kdiv 6 Chr ðkukrþ1 þ kpkrþ1 Þ:
¼ að2j að1 þ C 2pf ÞÞkrv k20 :
Applying the triangle inequality,
jku uh k1 þ jkp ph kdiv 6 jðku uI k1 þ kuI uh k1 Þ þ jðkp pI kdiv þ kpI ph kdiv Þ
j > 0, we obtain Choosing a ¼ 1þC 2 pf
j2 1 þ C 2pf
krv k20 ;
r
6 Ch ðkukrþ1 þ kpkrþ1 Þ: This completes the proof. h
which implies
kv k21 ¼ kv k20 þ krv k20 6 ð1 þ C 2pf Þkrv k20 6
which implies
a2 krv k20 P 2ajkrv k20 a2 C 2pf krv k20 a2 krv k20
Fððv ; qÞ; 0Þ P
185
ð1 þ C 2pf Þ2
j2
Fððv ; qÞ; 0Þ:
ð8Þ
On the other hand, we have
6 2ðFððv ; qÞ; 0Þ þ ð1 þ C 2pf ÞFððv ; qÞ; 0ÞÞ ¼ 2ð2 þ C 2pf ÞFððv ; qÞ; 0Þ;
ð9Þ
and there exists some constant C > 0 such that
kr qk20 6 2kr q þ a rv k20 þ 2ka rv k20 k20
6
C
j
2
Fððv ; qÞ; 0Þ:
ð10Þ
As a consequence of the coercivity estimate, one can easily show that problem (6) is uniquely solvable and the matrix A in the corresponding linear system An ¼ b is symmetric and positive definite. Theorem 2. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . Then the solution ðuh ; ph Þ of problem (6) satisfies the following estimate: there exists a constant C > 0 independent of j and h such that
jku uh k1 þ jkp ph kdiv 6 Chr ðkukrþ1 þ kpkrþ1 Þ:
ð11Þ
Proof. Let uI 2 Vh and pI 2 Qh be the standard finite element interpolants of u and p, respectively. Then from the interpolation theory, we have r
uð1Þ ¼ 0:
r
and kp pI kdiv 6 Ch kpkrþ1 ;
ð12Þ
for some constant C > 0 independent of j and h. Moreover, the error ðu uh ; p ph Þ has the following orthogonality property:
Bððuh u; ph pÞ; ðv h ; qh ÞÞ ¼ 0 8 ðv h ; qh Þ 2 Vh Qh :
ð13Þ
ð14Þ
In Fig. 1 the solution ðuh ; ph Þ of primitive LSFEM with piecewise linear finite elements ðr ¼ 1Þ is shown for j ¼ 102 and uniform mesh with h ¼ 1=32. Evidently, for such a relatively coarse mesh, the primitive LSFEM produces solution ðuh ; ph Þ that exhibit oscillatory behavior near the boundaries x ¼ 0; 1. In this case,
mesh Peclet number :¼
Now, combining (8)–(10) yields the conclusion. This completes the proof. h
ku uI k1 6 Ch kukrþ1
ju00 þ u0 ¼ 0 in ð0; 1Þ; uð0Þ ¼ 1;
kqk20 6 2ðkq þ jrv k20 þ j2 krv k20 Þ
6 2Fððv ; qÞ; 0Þ þ Ckrv
While error estimate (11) is optimal in the H10 ðXÞ Hðdiv; XÞnorm, the primitive LSFEM (6) does not yield satisfactory results if j is sufficiently small. In order to illustrate this fact, we consider the following simple one-dimensional problem:
jajh ¼ 1:5625 > 1: 2j
ð15Þ
That is, the primitive least-squares finite element approximation for this situation is convection-dominated. 3. The stabilized LSFEM In order to smear out the peaks occurring near the boundaries, in this section we will introduce the so-called stabilized LSFEM. As in the Galerkin least-squares method, we simply add a stabilized term in the previous L2 least-squares energy functional. However, due to the low regularity of functions in the space H10 ðXÞ Hðdiv; XÞ, we restrict the modified least-squares energy functional defined on the finite element space Vh Qh . Define Fs : Vh Qh ! R by
Fs ððv h ; qh Þ; f Þ ¼ kr qh þ a rv h f k20 þ kqh þ jrv h 0k20 X þ sT k jDv h þ a rv h f k20;T ; T2Th
where k k0;T denotes the L2 -norm restricted on the element T and sT > 0 is a parameter that need to be chosen appropriately. Taking the variation as in Section 2, we have the stabilized LSFEM:
Determine ðuh ; ph Þ 2 Vh Qh such that Bs ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Ls ððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Qh ;
ð16Þ
186
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p LSFEM ph
1
1
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u LSFEM u
0.2
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
x
0.6
0.8
1
x
Fig. 1. Exact solution vs. primitive least-squares finite element solution: h ¼ 1=32; j ¼ 102 .
where the bilinear form Bs and linear form Ls are defined as follows:
Hence,
Fs ððv h ; qh Þ; 0Þ
Bs ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 þ ðph þ jruh ; qh þ jrv h Þ0 X þ sT ðjDuh þ a ruh ; jDv h þ a rv h Þ0;T ;
P
X
sT ðf ; jDv h þ a rv h Þ0;T :
T2Th
Clearly, the stabilized LSFEM is a consistent formulation since (16) is satisfied with ðuh ; ph Þ replaced by the exact solution ðu; pÞ of problem (2). For deriving the coercivity estimate of Bs on the finite element space Vh Qh , we assume that the inverse inequality holds: there e > 0 independent of h such that exists a constant C
X
e C
2
hT kDuh k20;T 6 kruh k20
8 uh 2 Vh :
ð17Þ
T2Th
The inverse estimate (17) holds if the regular finite element partition Th of the domain X is quasi-uniform [25], namely, there exists b > 0 independent of h such that a constant C
b T h 6 Ch
e 2 Ch T 2ð1 þ
C 2pf Þ
C3
j2 kv
2 h k1
> 0;
ð19Þ
2 kqh k2div
þj
2
þ h ka rv
2 h k0
6 Fs ððv h ; qh Þ; 0Þ;
ð20Þ
Fs ððv h ; qh Þ; 0Þ P
T2Th
j2 2ð1 þ C 2pf Þ
2
krv h k20 þ Ch ka rv h k20 :
The remaining part of the proof is similar to that for Theorem 1.
h
Theorem 4. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameter sT is defined as in (19), then the solution ðuh ; ph Þ of problem (16) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that r
6 Ch ðkukrþ1 þ kpkrþ1 Þ:
Remark 1. Assume that Vh consists of piecewise linear polynomials ðr ¼ 1Þ. Taking sT ¼ s for any s > 0 and for all T 2 Th , we can derive the following estimate without assuming the inverse estimate (17):
C3
j2 kv h k21 þ j2 kqh k2div þ ska rv h k20 6 Fs ððv h ; qh Þ; 0Þ; pffiffiffi
Proof. Following the proof of Theorem 1, we can show that
1 þ C 2pf
j2 1þ
C 2pf
ð22Þ
based on which we have the following similar error estimate:
r1
j2
ð21Þ
Proof. The proof is similar to that for Theorem 2 with some additional modifications. We omit the details. h
6 Ch
¼
sT k jDv h k20;T 2k jDv h k20;T
jku uh k1 þ jkp ph kdiv þ ska rðu uh Þk0
for all ðv h ; qh Þ 2 Vh Qh .
Fs ððv h ; qh Þ; 0Þ P
X
Applying (17)–(19), we obtain
ð18Þ
then we have the following coercivity estimate: there exists a constant C 3 > 0 independent of j and h such that
krv h k20 þ
jku uh k1 þ jkp ph kdiv þ hka rðu uh Þk0
8 T 2 Th :
Theorem 3. Assume that the assumption (18) holds. If for each element T in Th we take
sT ¼
1þ
C 2pf
1 ka rv h k20;T þ ka rv h k20;T 2 X j2 1 2 2 2 2 ¼ r v þ s j k D v k þ r v k ka k k T h 0 h 0;T h 0;T : 2 1 þ C 2pf T2Th
T2Th
Ls ððv h ; qh ÞÞ ¼ ðf ; r qh þ a rv h Þ0 þ
j2
krv h k20 þ
X
sT k jDv h þ a rv h k20;T
T2Th
krv h k20 þ
X
sT k jDv h k20;T
T2Th
þ 2ðjDv h ; a rv h Þ0;T þ ka rv h k20;T :
pffiffiffi pffiffiffi ððh þ j s þ h sÞkukrþ1 þ hkpkrþ1 Þ:
ð23Þ
We now give the numerical simulation for the simple one-dimensional problem (10) with the same setting as that used for the primitive LSFEM in Fig. 1. The results of the stabilized LSFEM with sT ¼ s :¼ 0:085 for all T 2 Th are depicted in Fig. 2. We observe that the boundary peaks seem to disappear, but the height of uh is lower than the primitive least-squares finite element solution. Motivated by this interesting observation, we perform the numerical simulation with sT ¼ s :¼ 0:085 for all T 2 Th . The approximate solution we obtain is presented in Fig. 3. Note that
187
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p stabilized LSFEM ph
1
1
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u stabilized LSFEM u
0.2
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 2. Exact solution vs. stabilized least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; s ¼ 0:085.
the peaks still appear near the boundaries, but away from the boundary the height of uh coincides with the exact solution very well. This observation is a motivation of the heuristic scheme called the negatively stabilized streamline diffusion LSFEM developed in Section 5.
Theorem 5. Assume that the assumption (18) holds. If for each element T in Th we take
dT ¼
e 2 Ch T
j
< 0;
then we have the following coercivity estimate: there exists a constant C 4 > 0 independent of j and h such that
4. The streamline diffusion LSFEM
j kv 2
2 h k1
j2 h2 2 þ kqh k2div þ h ka rv h k20 j2 þ j2 h2 þ h4
There is another way to smear out the peaks occurring near the boundaries for the primitive least-squares finite element solution. The underlying idea is similar to that for the streamline diffusion Galerkin method (see, e.g., [25]). We introduce another first-order system formulation for problem (1) which induces the so-called streamline diffusion (SD) least-squares finite element method [22,23]. Introducing the new dependent variable,
C4
p :¼ jru þ daðjDu þ a ruÞ;
Fd ððv h ; qh Þ; 0Þ
ð24Þ
Fd ððv h ; qh Þ; f Þ 2
¼ kr qh þ a rv h fd k0 þ
X
2
kqh þ jrv h dT aðjDv h þ a rv h Þk0;T ;
T2Th
Proof. Let a > 0 be a constant that will be determined later. For any ðv h ; qh Þ 2 Vh Qh , we have ¼ kr qh þ a rv h av h k20 þ 2aðr qh þ a rv h ; v h Þ0 a2 kv h k20 X þ kqh þ jrv h dT aðjDv h þ a rv h Þ arv h k20;T T2Th
X
þ
2aðqh þ jrv h dT aðjDv h þ a rv h Þ; rv h Þ0;T a2 krv h k20 :
T2Th
Employing the Green formula (3) and noting that ða rv h ; v h Þ0 ¼ 0, we obtain
Fd ððv h ; qh Þ; 0Þ P a2 kv h k20 þ 2ajkrv h k20
where the small parameter d which introduce a correction in the streamline direction is chosen element-wise and denoted by dT . Taking the variation, we obtain the so-called streamline diffusion LSFEM:
X
X
2adT ðjDv h ; a rv h Þ0;T
T2Th
2adT ka rv h k20;T a2 krv h k20 :
T2Th
Since for any positive constant c,
Determine ðuh ; ph Þ 2 Vh Qh such that Bd ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Ld ððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Vh ;
ð28Þ
for all ðv h ; qh Þ 2 Vh Qh .
ð25Þ
We define the new least-squares energy functional Fd : Vh Qh ! R by
!
6 Fd ððv h ; qh Þ; 0Þ;
where d < 0 is an appropriate streamline diffusion parameter need to be determined later, then we have
r p þ a ru ¼ f þ dr ðf aÞ :¼ fd :
ð27Þ
ð26Þ
where the bilinear form Bd and linear form Ld are defined as follows, respectively:
Bd ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 X þ ðph þ jruh dT aðjDuh þ a ruh Þ; T2Th
X
P
Note that the streamline diffusion LSFEM is still a consistent scheme. For the streamline diffusion LSFEM (26), we have the following coercivity estimate:
X
1
adT j ck Dv h k20;T þ ka rv h k20;T ; c T2Th
we get that
Fd ððv h ; qh Þ; 0Þ
X
P 2ajkrv h k20 þ
qh þ jrv h dT aðjDv h þ a rv h ÞÞ0;T Ld ððv h ; qh ÞÞ ¼ ðfd ; r qh þ a rv h Þ0 :
2adT ðjDv h ; a rv h Þ0;T
T2Th
X T2Th
adT jckDv h k20;T
T2Th
adT 2
j ka rv h k20;T a2 ð1 þ C 2pf Þkrv h k20 ; c
where the Poincaré–Friedrichs inequality (4) has been used. Now j choosing a ¼ 2ð1þC 2 ; c ¼ j and using (27) with inverse estimate Þ pf (17), we obtain
188
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2 exact u negatively stabilized LSFEM u
0.2
exact p negatively stabilized LSFEM p
h
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 3. Exact solution vs. negatively stabilized least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; s ¼ 0:085.
Fd ððv h ; qh Þ; 0Þ P
j2 4ð1 þ C 2pf Þ
krv h k20 þ
X
e 2 Ch T
T2Th
2ð1 þ C 2pf Þ
ka rv h k20;T : ð29Þ
On the other hand,
kqh k20 6 2
X
X
kjrv h dT aðjDv h þ a rv h Þk20;T
T2Th
6 2Fd ððv h ; qh Þ; 0Þ þ 4j2 krv h k20 þ 8 þ8
X
j2 kv h k21 þ
jjdj jjdj þ 2d2 þ 1
! kqh k2div þ jjdjka rv h k20
6 Fs ððv h ; qh Þ; 0Þ;
ð33Þ
which implies the following error estimate:
kqh þ jrv h dT aðjDv h þ a rv h Þk20;T
T2Th
þ2
C4
X
d2T jaj2 j2 kDv h k20;T
T2Th
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi jjdj jku uh k1 þ kp ph kdiv þ jjdjka rðu uh Þk0 2 jjdj þ 2d þ 1 pffiffiffiffiffiffiffiffiffi r1 ð34Þ 6 Ch ðjjdj þ h jjdj þ hjdj þ hÞkukrþ1 þ hkpkrþ1 : In Fig. 4 we report the approximate solution ðuh ; ph Þ of the onedimensional problem (10) corresponding to the streamline diffusion LSFEM (23) with dT ¼ d :¼ 0:007 for all T 2 Th and the same values of j and h as before. Similar to the results obtained by stabilized LSFEM (16), boundary peaks disappear.
d2T jaj2 ka rv h k20;T :
T2Th
From (27) with inverse estimate (17), 2
kqh k20 6 2Fd ððv h ; qh Þ; 0Þ þ 4j2 krv h k20 þ Ch krv h k20 ! 2 2 h h þ C 2 Fd ððv h ; qh Þ; 0Þ ¼ C 1 þ 2 Fd ððv h ; qh Þ; 0Þ;
j
5. The negatively stabilized streamline diffusion LSFEM
ð30Þ
j
for some C > 0. Moreover, by the triangle inequality,
kr qh k20 6 2Fd ððv h ; qh Þ; 0Þ þ 2ka rv h k20 1 6 C 1 þ 2 Fd ððv h ; qh Þ; 0Þ: h Finally, the desired inequality (28) follows from (29)–(31).
ð31Þ
In the numerical results of streamline diffusion LSFEM reported in previous section, we find that the approximate solution is smooth but the height of uh is lower than the exact solution u. Motivated by the observation of the numerical results of the negatively stabilized LSFEM reported in Section 3, we now consider the streamline diffusion LSFEM with a negative stabilized term. Define the least-squares functional Fds : Vh Qh ! R by
Fds ððv h ; qh Þ; f Þ ¼ kr qh þ a rv h fd k20 þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2 h2 kp ph kdiv þ hka rðu uh Þk0 jku uh k1 þ j2 þ j2 h2 þ h4 ! 2 h r ð32Þ 6 Ch ð1 þ Þkukrþ1 þ kpkrþ1 :
j
Proof. The proof is similar to that for Theorem 2.h Remark 2. As in Remark 1, if Vh consists of piecewise linear polynomials ðr ¼ 1Þ and we take dT ¼ d for any d < 0 and for all T 2 Th , then the proof of Theorem 5 can be simplified to get the following estimate without assuming the inverse inequality (17):
kqh þ jrv h
T2Th
h
Theorem 6. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameter dT is defined as in (27), then the solution ðuh ; ph Þ of problem (26) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that
X
dT aðjDv h þ a rv h Þk20;T X þ sT k jDv h þ a rv h f k20;T ; T2Th
where dT < 0 and sT < 0 for all T 2 Th . Taking the variation, we get a negatively stabilized streamline diffusion LSFEM:
Determine ðuh ; ph Þ 2 Vh Qh such that Bds ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Lds ððv h ; qh ÞÞ 8ðv h ; qh Þ 2 Vh Vh ;
ð35Þ
where the bilinear form Bds and linear form Lds are defined as follows, respectively:
Bds ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 X þ ðph þ jruh dT aðjDuh þ a ruh Þ; qh T2Th
þ jrv h dT aðjDv h þ a rv h ÞÞ0;T X þ sT ðjDuh þ a ruh ; jDv h þ a rv h Þ0;T ; T2Th
189
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p SD LSFEM p
1
1
h
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u SD LSFEM uh
0.2 0
0
0.2
0
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 4. Exact solution vs. streamline diffusion least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; d ¼ 0:007.
Lds ððv h ; qh ÞÞ ¼ ðfd ; r qh þ a rv h Þ0 þ
X
sT ðf ; jDv h þ a rv h Þ0;T :
Proof. The proof is similar to that for Theorem 6. h
T2Th
Still, the negatively stabilized streamline diffusion LSFEM (35) is a consistent scheme and we have the following estimate: Theorem 7. Assume that the assumption (17) holds and for each element T in Th the parameter dT is defined as in (27). If take
sT ¼
e 2 Ch T eÞ 8ð1 þ C 2pf Þð1 þ C
< 0;
ð36Þ
then we have the following coercivity estimate: there exists a constant C 5 > 0 independent of j and h such that
C5
j2 kv h k21 þ
j2 h2 2 kqh k2div þ h ka rv h k20 2 j þ j2 h2 þ h4
!
6 Fds ððv h ; qh Þ; 0Þ;
ð37Þ
C5
j2 kv h k21 þ
!
jjdj jjdj þ 2d2 þ 1
kqh k2div þ jjdjka rv h k20
6 Fds ððv h ; qh Þ; 0Þ:
ð40Þ
Based on (40), we can derive the error estimate: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi jjdj kp ph kdiv þ jjdjka rðu uh Þk0 jku uh k1 þ 2 jjdj þ 2d þ 1 qffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi r1 6 Ch ð jjdj þ jdj þ 1Þh þ jjdj þ j3 jdj kukrþ1 þ hkpkrþ1 :
ð41Þ
for all ðv h ; qh Þ 2 Vh Qh . Proof. Similar to the proofs of Theorems 3 and 5, one can derive that
Fds ððv h ; qh Þ; 0Þ P
Remark 3. Assume that Vh consists of piecewise linear polynomials ðr ¼ 1Þ. If for each element T in Th we take dT ¼ d for any d < 0 jd < 0, then we have the following estimate withand sT ¼ s :¼ 1þC 2 pf out assuming the inverse inequality (17):
j2 C 2pf Þ
4ð1 þ X þ2
krv h k20 þ
X
e h2 C T
T2Th
2ð1 þ C 2pf Þ
ka rv h k20;T
sT ðk jDv h k20;T þ ka rv h k20;T Þ:
In Fig. 5 it is shown the negatively stabilized streamline diffusion least-squares finite element solution with dT ¼ d :¼ 0:007 and sT ¼ s :¼ 0:085 for all T 2 Th . Although the magnitude of s is larger than our estimate given in Remark 3, but it seems to achieve the optimal result for mesh size h ¼ 1=32. 6. Numerical experiments for 2D problems
T2Th
From (36), we conclude that
Fds ððv h ; qh Þ; 0Þ P
j2 8ð1 þ
C 2pf Þ
krv h k20 þ
X
e h2 C T
T2Th
4ð1 þ C 2pf Þ
ka rv h k20;T :
In this section, we perform a series of numerical simulations for four test problems by using the developed methods, including the primitive LSFEM (6), stabilized LSFEM (16), streamline diffusion LSFEM (26), and the negatively stabilized streamline diffusion LSFEM (35).
ð38Þ Moreover, we still have (30) and (31) with Fd ððv h ; qh Þ; 0Þ replaced by Fds ððv h ; qh Þ; 0Þ, and then the desired estimate (37) follows from (38), (30) and (31). h Similarly, we have the following error estimate: Theorem 8. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameters dT is defined as in (27) and sT as in (36), then the solution ðuh ; ph Þ of problem (35) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2 h2 jku uh k1 þ kp ph kdiv þ hka rðu uh Þk0 j2!þ j2 h2 þ h4 ! 2 h r kukrþ1 þ kpkrþ1 : ð39Þ 6 Ch 1þ
j
6.1. An example with analytic solution We first consider a simple problem with a ¼ ð0; 1Þ> on the unit square X ¼ ð0; 1Þ ð0; 1Þ subjected to the following boundary conditions [8]:
8 > < x ¼ 0 and 0 6 y 6 1; u ¼ 0 when x ¼ 1 and 0 6 y 6 1; > : y ¼ 1 and 0 6 x 6 1;
u ¼ sinðpxÞ when y ¼ 0 and 0 6 x 6 1: One can check that the exact solution is given by
u¼
1 ðem2 m1 em1 y em2 y Þ sinðpxÞ; em2 m1 1
where
m1 ¼
1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4j2 p2 2j
and m2 ¼
1þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4j2 p2 : 2j
190
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2 exact u negatively stabilized SD LSFEM uh
0.2 0
0
0.2
0.4
0.6
0.8
exact p negatively stabilized SD LSFEM p
h
0 1
−0.2
0
0.2
x
0.4
0.6
0.8
1
x
Fig. 5. Exact solution vs. negatively stabilized streamline diffusion least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; d ¼ 0:007; s ¼ 0:085.
Table 1 Relative error behavior for
j ¼ 1 and r ¼ 1.
Method
u=p
h ¼ 1=32
h ¼ 1=64
h ¼ 1=128
h ¼ 1=256
LSFEM
u
2
L
9.3538E4
2.3849E4
6.0488E5
1.5344E5
1.98
u
H1
2.0763E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9168E2
7.2544E3
2.7362E3
1.0290E3
1.41
L2
4.9617E2 9.1326E4
2.5847E2 2.3295E4
1.3219E2 5.9100E5
6.6925E3 1.5001E5
0.96 1.98
u
H1
2.0762E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9167E2
7.2544E3
2.7362E3
1.0290E3
1.41
L2
4.9617E2 9.2450E4
2.5847E2 2.3577E4
1.3219E2 5.9807E5
6.6925E3 1.5176E5
0.96 1.98
u
H1
2.0762E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9166E2
7.2542E3
2.7361E3
1.0290E3
1.41
L2
4.9617E2 9.3343E4
2.5847E2 2.3801E4
1.3219E2 6.0372E5
6.6925E3 1.5314E5
0.96 1.98
p u
Stabilized LSFEM
p u
SD LSFEM
p u
Negatively stabilized
Norm
Order
u
H1
2.0763E2
1.0379E2
5.1895E3
2.5947E3
1.00
SD
p
1.9166E2
7.2541E3
2.7360E3
1.0290E3
1.41
LSFEM
p
L2 HðdivÞ
4.9617E2
2.5847E2
1.3219E2
6.6925E3
0.96
(Stabilized LSFEM:
2
2
2
2
s ¼ 0:1h ; SD LSFEM: d ¼ 0:05h =j; negatively stabilized SD LSFEM: d ¼ 0:07h =j; s ¼ 0:06h ).
Table 2 Relative error behavior for
j ¼ 102 and r ¼ 1.
Method
u=p
Norm
h ¼ 1=32
h ¼ 1=64
h ¼ 1=128
h ¼ 1=256
Order
LSFEM
u
L2
9.7176E2
3.0693E2
8.2541E3
2.1150E3
1.84
u
H1
4.4607E5
2.9571E5
1.6109E5
8.2399E2
0.81
p
L2 HðdivÞ
3.8496E5
1.2947E5
3.6118E2
9.7787E3
1.77
L2
4.9399E5 1.7826E5
3.2868E5 7.0424E2
1.7922E5 2.0870E2
9.1710E2 5.4864E3
0.81 1.67
u
H1
4.3500E5
2.9518E5
1.6126E5
8.2437E2
0.80
p
L2 HðdivÞ
3.9890E5
1.3585E5
3.8033E2
1.0299E2
1.76
L2
4.7611E5 6.8719E2
3.2665E5 2.7989E2
1.7917E5 8.3983E3
9.1720E2 2.2232E3
0.79 1.65
u
H1
4.1448E5
2.8767E5
1.5982E5
8.2231E2
0.78
p
L2 HðdivÞ
5.1168E5
1.6527E5
4.5548E2
1.2106E2
1.80
L2
4.5991E5 4.3060E2
3.1980E5 1.3619E2
1.7781E5 3.6980E3
9.1522E2 9.6171E4
0.78 1.83
u
H1
4.1941E5
2.8944E5
1.6012E5
8.2270E2
0.78
p
L2 HðdivÞ
5.8724E5
1.8574E5
5.0631E2
1.3325E2
1.82
4.6612E5
3.2198E5
1.7818E5
9.1571E2
0.78
p u
Stabilized LSFEM
p u
SD LSFEM
Negatively stabilized SD LSFEM (Stabilized LSFEM:
p u
p
s ¼ 0:1h2 ; SD LSFEM: d ¼ 0:05h2 =j; negatively stabilized SD LSFEM: d ¼ 0:07h2 =j; s ¼ 0:06h2 ).
Substituting the above exact solution into problem (2), we can readily get the right-hand side function f. We consider two different values for the diffusivity, namely, j ¼ 1 and j ¼ 102 . For simplifying the numerical implementation, we shall assume that the unit square domain X is uniformly 2 partitioned into a set of h square subdomains with side-length
h ¼ 2‘ for ‘ ¼ 5; 6; 7; 8. We will apply continuous piecewise bilinear finite elements ðr ¼ 1Þ to approximate all components of the exact solution for all least-squares finite element schemes. In all computations, a double precision conjugate gradient solver is applied to solve the associated linear systems. We compute the relative error in the L2 ðXÞ-norm and H1 ðXÞ-norm for u, and L2 ðXÞ-norm
191
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact solution
exact solution 1 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0
y
0.5
x
0
1 1
0
0.2
Fig. 6. Exact solution u and its contours:
LSFEM
0.4
0.6
0.8
1
j ¼ 102 .
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
1 1
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 7. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y 1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 8. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; s ¼ 0:1h .
SD LSFEM
SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 9. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; d ¼ 0:05h =j.
and Hðdiv; XÞ-norm for p. The results are collected in Tables 1 and 2. The elevation and contour plots for the exact solution u and
approximate solutions uh for h ¼ 1=32 and in Figs. 6–10.
j ¼ 102 are displayed
192
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y 1 1
0.5
x
0
0
0.2
0.4
0.6
0.8
1 2
2
Fig. 10. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; d ¼ 0:07h =j; s ¼ 0:06h .
LSFEM
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0
0.5 y 1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 6
Fig. 11. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x 0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 12. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 106 ; s ¼ 0:2h .
SD LSFEM
SD LSFEM
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x 0.5
0
0
0.2
0.4
0.6
0.8
1 6
2
Fig. 13. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 ; d ¼ 0:00002h =j.
Numerical results reported in Table 1 for j ¼ 1 confirm the theoretical error estimates with r ¼ 1. That is, all the primitive LSFEM and its variants developed in this paper achieve optimal conver-
gence in the H1 ðXÞ-norm for u and in the Hðdiv; XÞ-norm for p. Furthermore, an examination of the numerical results given in Table 1 shows that the convergence rate for u in L2 ðXÞ-norm seems to be
193
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
1 1
0.5
0
0
0.2
0.4
0.6
0.8
1 2
2
Fig. 14. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 106 ; d ¼ 0:0000163h =j; s ¼ 0:06h .
LSFEM
LSFEM 1 0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
1 0.8 0.6 0.4 0.2 0
0 0 0.5
y
0.5
x
0
1 1
0
0.2
0.4
0.6
0.8
1 4
Fig. 15. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
1 0.8 0.6 0.4 0.2 0
0 0 0.5
y
0.5
x
0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 16. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; s ¼ 1.
SD LSFEM
SD LSFEM 1
1 0.8
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.6 0.4 0.2 0
0 0 0.5
0.5 y
1 1
x
0
0 0
0.2
0.4
0.6
0.8
1
Fig. 17. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 104 .
optimal, although we have not provided the analysis. Also, from Table 2 for j ¼ 102 , it is no surprise that the convergence behavior
for all unknowns is getting worse when the diffusivity j is getting small.
194
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1 0.8
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.6 0.4 0.2 0
0 0 0.5
y
x
0.5 0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 18. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 104 ; s ¼ 0:5.
LSFEM
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
0.5
0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 19. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
0.5
1 1
0
0 0
0.2
0.4
0.6
0.8
1 2
Fig. 20. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; s ¼ 0:2h .
6.2. Boundary layer problem
6.3. Interior layer problem
We consider the test case where X ¼ ð0; 1Þ ð0; 1Þ; j ¼ 106 ; a ¼ ð0; 1Þ> , f ¼ 0, and the boundary conditions are given by
ary data with a uniform velocity field a ¼
u ¼ 1 on fðx; yÞ : 0 6 x 6 1; y ¼ 0g; u ¼ 0 on fðx; yÞ : 0 6 x 6 1; y ¼ 1g; and the homogeneous Neumann boundary condition are imposed on the other two sides, namely,
p1 ¼ 0 on fðx; yÞ : x ¼ 0; 1 and 0 6 y 6 1g: A uniform mesh with h ¼ 25 and piecewise bilinear finite elements are employed for all least-squares methods. The results are shown in Figs. 11–14. One can observe that the negatively stabilized streamline diffusion least-squares method is able to better capture the boundary layer behavior when compared with other leastsquares methods.
Let us consider the propagation of a discontinuity in the bound > p1ffiffi ; p1ffiffi on the unit 2 2
square domain X ¼ ð0; 1Þ ð0; 1Þ with f ¼ 0 and boundary conditions are given by
j ¼ 104 . The
u ¼ 1 on fðx; yÞ : x ¼ 0; 0 6 y 6 1g [ fðx; yÞ : 0 6 x 6 1; y ¼ 1g; u ¼ 0 on fðx; yÞ : 0 < x 6 1; y ¼ 0g [ fðx; yÞ : x ¼ 1; 0 6 y < 1g: We have an interior layer due to the discontinuity of the boundary data. A uniform mesh with h ¼ 25 and piecewise bilinear finite elements are employed for all least-squares methods. In Figs. 15–18, the elevation and contour plots for approximate solution uh produced by various methods are displayed. Except the primitive LSFEM, all other least-squares methods perform very well.
195
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
SD LSFEM
SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.8 0.6 0.4
0.4
0.4
0.2 0 0
0.2
0.2 0.5
0 y
x
0.5
0
1 1
0 0
0.2
0.4
0.6
0.8
1 2
Fig. 21. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 0:0017h =j.
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1 1
1
0.8 0.8
0.8 0.6 0.4
0.6
0.6
0.4
0.4
0.2
0.2
0.2 0 0 0.5
0 y 1 1
x
0.5
0
0 0
0.2
0.4
0.6
0.8
1 2
2
Fig. 22. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 0:0017h =j; s ¼ 0:055h .
6.4. Boundary and interior layer problem We now construct a model problem possessing both boundary and interior layer behavior on the unit square domain X ¼ ð0; 1Þ ð0; 1Þ. As in Section 6.3, we consider the propagation of a discontinuity in the boundary data with constant velocity field pffiffi> a ¼ 12 ; 23 of size one forming a 60° with the x-axis. In this example, f ¼ 0; j ¼ 104 , and the boundary conditions are given by
the primitive LSFEM performs poorly for the convection-dominated model while the other three variants perform considerably better for interior layer problems, and the negatively stabilized streamline diffusion LSFEM is able to better capture the boundary layer behavior when compared with other LSFEMs. But for problem possessing both interior and boundary layer structures in the solution, more investigation is needed in order to reach an efficient LSFEM. The study in this direction is in progress.
u ¼ 1 on fðx; yÞ : x ¼ 0; 0 6 y < 1g [ fðx; yÞ : 0 6 x 6 1=2; y ¼ 0g; u ¼ 0 on fðx; yÞ : 0 6 x 6 1; y ¼ 1g [ fðx; yÞ : x ¼ 1; 0 6 y 6 1g
Acknowledgment
[ fðx; yÞ : 1=2 < x 6 1; y ¼ 0g:
This work was partially supported by the National Science Council of Taiwan.
The elevation and contour plots of all approximate solutions are displayed in Figs. 19–22. Unfortunately, all the primitive least-squares method and its variants do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. 7. Concluding remarks In this paper, we have proposed and analyzed the primitive LSFEM and its three variants, including the stabilized LSFEM, streamline diffusion LSFEM and the negatively stabilized streamline diffusion LSFEM. The coercivity (stability) estimates for all methods are given. One of the common advantages of these least-squares finite element methods is that the resulting linear system is symmetric and positive definite. It was observed that
References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convective dominated flows with a particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Engrg. 32 (1982) 199–259. [2] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective–diffusive model, Comput. Meth. Appl. Mech. Engrg. 95 (1992) 253–276. [3] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Meth. Appl. Mech. Engrg. 73 (1989) 173–189. [4] L.P. Franca, G. Hauke, A. Masud, Revisiting stabilized finite element methods for the advective–diffusive equation, Comput. Meth. Appl. Mech. Engrg. 195 (2006) 1560–1572. [5] F. Brezzi, A. Russo, Choosing bubbles for advection–diffusion problems, Math. Mod. Meth. Appl. Sci. 4 (1994) 571–587.
196
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
[6] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Meth. Appl. Mech. Engrg. 166 (1998) 25–33. [7] F. Brezzi, T.J.R. Hughes, L.D. Marini, A. Russo, E. Süli, A priori error analysis of residual-free bubbles for advection–diffusion problems, SIAM J. Numer. Anal. 36 (1999) 1933–1948. [8] L.P. Franca, A. Nesliturk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by a two-level finite element method, Comput. Meth. Appl. Mech. Engrg. 166 (1998) 35–49. [9] P.B. Bochev, M.D. Gunzburger, Finite element methods of least-squares type, SIAM Rev. 40 (1998) 789–837. [10] B.-N. Jiang, The Least-Squares Finite Element Method, Springer-Verlag, Berlin, 1998. [11] P.B. Bochev, M.D. Gunzburger, Analysis of least-squares finite element methods for the Stokes equations, Math. Comp. 63 (1994) 479–506. [12] Z. Cai, R. Lazarov, T.A. Manteuffel, S.F. McCormick, First-order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal. 31 (1994) 1785–1799. [13] Z. Cai, T.A. Manteuffel, S.F. McCormick, First-order system least squares for second-order partial differential equations: Part II, SIAM J. Numer. Anal. 34 (1997) 425–454. [14] C.L. Chang, B.-N. Jiang, An error analysis of least-squares finite element method of velocity–pressure–vorticity formulation for Stokes problem, Comput. Meth. Appl. Mech. Engrg. 84 (1990) 247–255. [15] C.L. Chang, S.-Y. Yang, C.-H. Hsu, A least-squares finite element method for incompressible flow in stress–velocity–pressure version, Comput. Meth. Appl. Mech. Engrg. 128 (1995) 1–9.
[16] C.L. Chang, S.-Y. Yang, Analysis of the L2 least-squares finite element method for the velocity–vorticity–pressure Stokes equations with velocity boundary conditions, Appl. Math. Comput. 130 (2002) 121–144. [17] J.M. Deang, M.D. Gunzburger, Issues related to least-squares finite element methods for the Stokes equations, SIAM J. Sci. Comput. 20 (1998) 878–906. [18] A.I. Pehlivanov, G.F. Carey, R.D. Lazarov, Least-squares mixed finite elements for second-order elliptic problems, SIAM J. Numer. Anal. 31 (1994) 1368–1377. [19] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994. [20] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, SpringerVerlag, New York, 1991. [21] J.M. Fiard, T.A. Manteuffel, S.F. McCormick, First-order system least squares (FOSLS) for convection–diffusion problems: numerical results, SIAM J. Sci. Comput. 19 (1998) 1958–1979. [22] R.D. Lazarov, L. Tobiska, P.S. Vassilevski, Streamline diffusion least-squares mixed finite element methods for convection–diffusion problems, East-West J. Numer. Math. 5 (1997) 249–264. [23] R.D. Lazarov, P.S. Vassilevski, Least-squares streamline diffusion finite element approximations to singularly perturbed convection–diffusion problems, in: L.G. Vulkov, J.J.H. Miller, G.I. Shishkin (Eds.), Analytical and Numerical Methods for Singularly Perturbed Problems, Nova Science Publishing House, 2000, pp. 83–94. [24] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, New York, 1986. [25] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, Cambridge, 1987.
Contents lists available at ScienceDirect
Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma
On efficient least-squares finite element methods for convection-dominated problems Po-Wen Hsieh, Suh-Yuh Yang * Department of Mathematics, National Central University, Jhongli City 32001, Taiwan
a r t i c l e
i n f o
Article history: Received 2 September 2008 Accepted 30 September 2009 Available online 9 October 2009 MSC: 65N12 65N15 65N30 Keywords: Convection-dominated problems Boundary and interior layers Finite element methods Least-squares Stabilized methods Streamline diffusion
a b s t r a c t This paper focuses on the least-squares finite element method and its three variants for obtaining efficient numerical solutions to convection-dominated convection–diffusion problems. The coercivity estimates for the corresponding homogeneous least-squares energy functionals are derived, and based on which error estimates are established. One of the common advantages of these least-squares methods is that the resulting linear system is symmetric and positive definite. Numerical experiments that demonstrate the theoretical analysis of the developed methods are presented. It was observed that the primitive least-squares method performs poorly for convection-dominated problems while the stabilized, streamline diffusion and negatively stabilized streamline diffusion least-squares methods perform considerably better for interior layer problems, and the negatively stabilized streamline diffusion leastsquares method is able to better capture the boundary layer behavior when compared with other least-squares methods. But all the least-squares methods do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction Let X be a bounded convex polygonal domain in R2 with boundary @ X. We consider the following convection–diffusion problem:
jDu þ a ru ¼ f u¼g
in X;
on @ X;
ð1Þ
where 0 < j 6 1 is the constant viscosity coefficient (or diffusivity), a ¼ ða1 ; a2 Þ> is the velocity field, f is a source function and g is a given boundary data. It is well-known that the solution u of problem (1) can exhibit localized phenomena such as boundary and interior layers, i.e., narrow regions where the derivative of u is very large. This is the case for instance if the problem is convection-dominated, i.e., 0 < j jaj. For simplifying the mathematical analysis, in what follows we assume that a is a constant vector with jaj ¼ 1 and g ¼ 0 on the boundary @ X. It is already known that conventional numerical methods are lacking in either stability or accuracy for the convection-dominated problems. For example, the standard Galerkin method using piecewise linear elements performs poorly for the convection-dominated
* Corresponding author. Tel.: +886 3 4227151x65130; fax: +886 3 4257379. E-mail address: [email protected] (S.-Y. Yang). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.09.029
model, in which large spurious oscillations are frequently exhibited. In order to defeat this difficulty, the subject of stabilized finite element methods has existed for more than twenty years and is still attractive today (see, e.g., [1–3] or recent review [4] by Franca et al.). The stabilized methods are formed by adding to the standard Galerkin method variational terms that are consistent and meshdependent. This additional term not only improves the numerical stability of the Galerkin method but also preserves good accuracy. Another possible alternative is still to apply the standard finite element method to solve the convection-dominated problem by enriching piecewise linear polynomials with bubble functions (see the early tutorial paper [5] by Brezzi and Russo, or [6–8] and many references therein). The bubble function is chosen so that the computed solution satisfies the original second-order differential equation in the interior of each element and vanishes on the boundary of each element. This produces a stabilized-like formulation. However, the resulting linear systems of above-mentioned two approaches are not symmetric. On the other hand, for the past two decades, least-squares finite element methods (LSFEMs) have become more and more frequently used to approximate the solution of first-order system of partial differential equations arising from fluid and solid mechanics. We refer the reader to [9,10] for some of the vast literature. In a broad sense, the least-squares approach is to minimize the
184
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
residual (measured in some norms) of a first-order system of partial differential equations over a suitable function space. In general, such a first-order system is obtained from an originally higher-order system by introducing one or more additional physically unknown variables such as vorticity, stress, or flux (see, e.g., [11–18]). By applying the least-squares principle on this extended first-order system, the smoothness requirement on the function space can then be relaxed. It is already known that the specific features of the least-squares finite element approach which give it advantages relative to, for instance, the mixed finite element approach [19,20], are as follows: it leads to a minimization problem; it is not subject to the Ladyzhenskaya–Babuska–Brezzi condition; simple equal low-order finite elements, such as the continuous linear elements, can be used for the approximation of all unknowns; the resulting linear system is symmetric and positive definite; the value of the homogeneous least-squares functional of the approximate solution provides a practical and sharp a posteriori error estimator at no additional cost, etc. Unfortunately, as we will see in this paper, the primitive LSFEM is not able to achieve good performance for convection-dominated problems. In [21], the authors proposed an exponentially weighted LSFEM to achieve H1 -like accuracy of the discretization and optimal multigrid convergence uniformly in the size of convection. However, their norms incorporate a scaling that has the effect of dramatically changing the influence of the boundary layer. The aim of this paper is to combine the advantages of abovementioned stabilized and least-squares finite element approaches to tackle the convection-dominated problems. We will introduce and analyze the primitive LSFEM and its three variants. To this aim, we have to recast problem (1) into the following first-order system problem by introducing the new dependent variable p ¼ jru in X:
product ðp; qÞdiv :¼ ðp; qÞ0 þ ðr p; r qÞ0 and norm kqk2div :¼ kqk20 þ kr qk20 . In the analysis of the developed methods, we will frequently use the following Green formula:
8 > < r p þ a ru ¼ f p þ jru ¼ 0 > : u¼0
Fððu; pÞ; f Þ ¼ 0 ¼ minfFððv ; qÞ; f Þ : ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞg:
in X; in X;
ð2Þ
on @ X:
We first apply the primitive LSFEM to approximate the solution to problem (2), in which the corresponding least-squares energy functional is defined to be the sum of the squared L2 norms of the residuals of the partial differential equations over an appropriate product space, and the least-squares finite element solution is defined to be the minimizer of the functional over a continuous piecewise polynomial finite element space. As we will see further below, the primitive LSFEM performs poorly for convection-dominated problems. We then introduce three variants of the primitive LSFEM that will be introduced in next several sections, including the stabilized LSFEM, the streamline diffusion LSFEM [22,23] and the negatively stabilized streamline diffusion LSFEM. Numerical simulations will be performed to demonstrate the theoretical analysis of the proposed methods. More importantly, it was observed that the stabilized, streamline diffusion and negatively stabilized streamline diffusion LSFEMs perform considerably better than the primitive LSFEM for interior layer problems, and the negatively stabilized streamline diffusion LSFEM is able to better capture the boundary layer behavior when compared with other least-squares methods. But regrettably, all the least-squares methods do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. Throughout this paper, we will use standard notation and definition for the Sobolev space Hm ðXÞ for non-negative integer m (cf. [19,24]). The standard associated inner product and norm are denoted by ð; Þm and k km , respectively. As usual, L2 ðXÞ ¼ H0 ðXÞ and H10 ðXÞ ¼ fv 2 H1 ðXÞ and v j@ X ¼ 0g. We also need the space Hðdiv; XÞ ¼ fq 2 ½L2 ðXÞ2 and r q 2 L2 ðXÞg with standard inner
ðq; rv Þ0 þ ðr q; v Þ0 ¼ ðq n; v Þ0;@ X
8 q 2 Hðdiv; XÞ and
v 2 H1 ðXÞ; ð3Þ
where n denotes the unit outward normal vector to @ X, and we also need the following Poincaré–Friedrichs inequality [24]: there exists a constant C pf > 0 such that
kv k0 6 C pf krv k0
8
v 2 H10 ðXÞ:
ð4Þ
The remainder of this paper is organized as follows. We introduce and analyze the primitive LSFEM, the stabilized LSFEM, the streamline diffusion LSFEM, and the negatively stabilized streamline diffusion LSFEM in Sections 2–5, respectively. The coercivity estimates for the corresponding homogeneous least-squares functionals are established, and based on which error estimates are derived. Numerical results for various test cases are reported in Section 6. Some concluding remarks are given in Section 7. 2. The primitive LSFEM Define the L2 least-squares energy functional F : H10 ðXÞ Hðdiv; XÞ ! R for the extended first-order problem (2) by
Fððv ; qÞ; f Þ ¼ kr q þ a rv f k20 þ kq þ jrv 0k20 : Note that the energy functional Fð; f Þ is defined to be the sum of the squared L2 -norms of the residuals on the product function space H10 ðXÞ Hðdiv; XÞ. Thus, if ðu; pÞ 2 H10 ðXÞ Hðdiv; XÞ is an exact solution of problem (2), then ðu; pÞ must be a zero minimizer of the functional Fð; f Þ on H10 ðXÞ Hðdiv; XÞ, namely,
Since Fððu; pÞ þ eðv ; qÞ; f Þ is a nonnegative quadratic functional in e 2 R for any given ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ, we have
d Fððu; pÞ þ eðv ; qÞ; f Þje¼0 ¼ 0; de which is equivalent to
Bððu; pÞ; ðv ; qÞÞ ¼ Lððv ; qÞÞ 8 ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ:
ð5Þ
Here the bilinear form Bð; Þ and the linear form LðÞ are defined as following, respectively:
Bððu; pÞ; ðv ; qÞÞ ¼ ðr p þ a ru; r q þ a rv Þ0 þ ðp þ jru; q þ jrv Þ0 ; Lððv ; qÞÞ ¼ ðf ; r q þ a rv Þ0 : Now we consider the finite element formulation of the leastsquares method (5) for problem (2). We introduce the finite element spaces Vh # H10 ðXÞ and Qh # Hðdiv; XÞ which consist of piecewise polynomials of degree less or equal than r (r P 1 integer) over a regular finite element partition Th of the domain X. We denote by hT the diameter of element T 2 Th and set h :¼ maxT2Th hT . The primitive LSFEM for problem (2) is then defined to be the following:
Determine ðuh ; ph Þ 2 Vh Qh such that Bððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Lððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Qh :
ð6Þ
Once the basis functions of the finite element space Vh Qh are chosen, problem (6) is equivalent to solve a linear system problem, An ¼ b, where unknown vector n consists the coordinates of the finite element solution ðuh ; ph Þ with respect to the chosen basis functions. The unique solvability of problem (6) and error estimates of the approximate solution ðuh ; ph Þ is mainly depending on the coercivity
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
(i.e., stability) and boundedness estimates of bilinear form B. Note that Bððv ; qÞ; ðv ; qÞÞ ¼ Fððv ; qÞ; 0Þ for all ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ. Theorem 1. There exist constants C 1 ; C 2 > 0 independent of such that
C1
j2 kv k21 þ j2 kqk2div
j and h
6 Fððv ; qÞ; 0Þ 6 C 2 kv k21 þ kqk2div ;
According to the orthogonality property (13), we have
Bððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ ¼ Bððuh u; ph pÞ þ ðu uI ; p pI Þ; ðuh uI ; ph pI ÞÞ ¼ Bððu uI ; p pI Þ; ðuh uI ; ph pI ÞÞ
ð7Þ
1
6 B2 ððu uI ; p pI Þ; ðu uI ; p pI ÞÞ 1
for all ðv ; qÞ 2 H10 ðXÞ Hðdiv; XÞ.
B2 ððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ;
Proof. The boundedness estimate can be easily proven. We proceed to show the validity of the coercivity estimate, namely, the left hand side of (7). Let a be a positive constant that will be determined later. Applying the Poincaré–Friedrichs inequality (4), we have
Fððv ; qÞ; 0Þ ¼ kr q þ a rv av k20 þ 2aðr q þ a rv ; v Þ0 a2 kv k20 þ kq þ jrv arv k20 þ 2aðq þ jrv ; rv Þ0
1
B2 ððuh uI ; ph pI Þ; ðuh uI ; ph pI ÞÞ 1
6 B2 ððu uI ; p pI Þ; ðu uI ; p pI ÞÞ: From above inequality, Theorem 1, interpolation property (12) and the boundedness of the bilinear form B on the space H10 ðXÞ Hðdiv; XÞ,
jkuh uI k1 þ jkph pI kdiv 6 Chr ðkukrþ1 þ kpkrþ1 Þ:
¼ að2j að1 þ C 2pf ÞÞkrv k20 :
Applying the triangle inequality,
jku uh k1 þ jkp ph kdiv 6 jðku uI k1 þ kuI uh k1 Þ þ jðkp pI kdiv þ kpI ph kdiv Þ
j > 0, we obtain Choosing a ¼ 1þC 2 pf
j2 1 þ C 2pf
krv k20 ;
r
6 Ch ðkukrþ1 þ kpkrþ1 Þ: This completes the proof. h
which implies
kv k21 ¼ kv k20 þ krv k20 6 ð1 þ C 2pf Þkrv k20 6
which implies
a2 krv k20 P 2ajkrv k20 a2 C 2pf krv k20 a2 krv k20
Fððv ; qÞ; 0Þ P
185
ð1 þ C 2pf Þ2
j2
Fððv ; qÞ; 0Þ:
ð8Þ
On the other hand, we have
6 2ðFððv ; qÞ; 0Þ þ ð1 þ C 2pf ÞFððv ; qÞ; 0ÞÞ ¼ 2ð2 þ C 2pf ÞFððv ; qÞ; 0Þ;
ð9Þ
and there exists some constant C > 0 such that
kr qk20 6 2kr q þ a rv k20 þ 2ka rv k20 k20
6
C
j
2
Fððv ; qÞ; 0Þ:
ð10Þ
As a consequence of the coercivity estimate, one can easily show that problem (6) is uniquely solvable and the matrix A in the corresponding linear system An ¼ b is symmetric and positive definite. Theorem 2. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . Then the solution ðuh ; ph Þ of problem (6) satisfies the following estimate: there exists a constant C > 0 independent of j and h such that
jku uh k1 þ jkp ph kdiv 6 Chr ðkukrþ1 þ kpkrþ1 Þ:
ð11Þ
Proof. Let uI 2 Vh and pI 2 Qh be the standard finite element interpolants of u and p, respectively. Then from the interpolation theory, we have r
uð1Þ ¼ 0:
r
and kp pI kdiv 6 Ch kpkrþ1 ;
ð12Þ
for some constant C > 0 independent of j and h. Moreover, the error ðu uh ; p ph Þ has the following orthogonality property:
Bððuh u; ph pÞ; ðv h ; qh ÞÞ ¼ 0 8 ðv h ; qh Þ 2 Vh Qh :
ð13Þ
ð14Þ
In Fig. 1 the solution ðuh ; ph Þ of primitive LSFEM with piecewise linear finite elements ðr ¼ 1Þ is shown for j ¼ 102 and uniform mesh with h ¼ 1=32. Evidently, for such a relatively coarse mesh, the primitive LSFEM produces solution ðuh ; ph Þ that exhibit oscillatory behavior near the boundaries x ¼ 0; 1. In this case,
mesh Peclet number :¼
Now, combining (8)–(10) yields the conclusion. This completes the proof. h
ku uI k1 6 Ch kukrþ1
ju00 þ u0 ¼ 0 in ð0; 1Þ; uð0Þ ¼ 1;
kqk20 6 2ðkq þ jrv k20 þ j2 krv k20 Þ
6 2Fððv ; qÞ; 0Þ þ Ckrv
While error estimate (11) is optimal in the H10 ðXÞ Hðdiv; XÞnorm, the primitive LSFEM (6) does not yield satisfactory results if j is sufficiently small. In order to illustrate this fact, we consider the following simple one-dimensional problem:
jajh ¼ 1:5625 > 1: 2j
ð15Þ
That is, the primitive least-squares finite element approximation for this situation is convection-dominated. 3. The stabilized LSFEM In order to smear out the peaks occurring near the boundaries, in this section we will introduce the so-called stabilized LSFEM. As in the Galerkin least-squares method, we simply add a stabilized term in the previous L2 least-squares energy functional. However, due to the low regularity of functions in the space H10 ðXÞ Hðdiv; XÞ, we restrict the modified least-squares energy functional defined on the finite element space Vh Qh . Define Fs : Vh Qh ! R by
Fs ððv h ; qh Þ; f Þ ¼ kr qh þ a rv h f k20 þ kqh þ jrv h 0k20 X þ sT k jDv h þ a rv h f k20;T ; T2Th
where k k0;T denotes the L2 -norm restricted on the element T and sT > 0 is a parameter that need to be chosen appropriately. Taking the variation as in Section 2, we have the stabilized LSFEM:
Determine ðuh ; ph Þ 2 Vh Qh such that Bs ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Ls ððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Qh ;
ð16Þ
186
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p LSFEM ph
1
1
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u LSFEM u
0.2
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
x
0.6
0.8
1
x
Fig. 1. Exact solution vs. primitive least-squares finite element solution: h ¼ 1=32; j ¼ 102 .
where the bilinear form Bs and linear form Ls are defined as follows:
Hence,
Fs ððv h ; qh Þ; 0Þ
Bs ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 þ ðph þ jruh ; qh þ jrv h Þ0 X þ sT ðjDuh þ a ruh ; jDv h þ a rv h Þ0;T ;
P
X
sT ðf ; jDv h þ a rv h Þ0;T :
T2Th
Clearly, the stabilized LSFEM is a consistent formulation since (16) is satisfied with ðuh ; ph Þ replaced by the exact solution ðu; pÞ of problem (2). For deriving the coercivity estimate of Bs on the finite element space Vh Qh , we assume that the inverse inequality holds: there e > 0 independent of h such that exists a constant C
X
e C
2
hT kDuh k20;T 6 kruh k20
8 uh 2 Vh :
ð17Þ
T2Th
The inverse estimate (17) holds if the regular finite element partition Th of the domain X is quasi-uniform [25], namely, there exists b > 0 independent of h such that a constant C
b T h 6 Ch
e 2 Ch T 2ð1 þ
C 2pf Þ
C3
j2 kv
2 h k1
> 0;
ð19Þ
2 kqh k2div
þj
2
þ h ka rv
2 h k0
6 Fs ððv h ; qh Þ; 0Þ;
ð20Þ
Fs ððv h ; qh Þ; 0Þ P
T2Th
j2 2ð1 þ C 2pf Þ
2
krv h k20 þ Ch ka rv h k20 :
The remaining part of the proof is similar to that for Theorem 1.
h
Theorem 4. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameter sT is defined as in (19), then the solution ðuh ; ph Þ of problem (16) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that r
6 Ch ðkukrþ1 þ kpkrþ1 Þ:
Remark 1. Assume that Vh consists of piecewise linear polynomials ðr ¼ 1Þ. Taking sT ¼ s for any s > 0 and for all T 2 Th , we can derive the following estimate without assuming the inverse estimate (17):
C3
j2 kv h k21 þ j2 kqh k2div þ ska rv h k20 6 Fs ððv h ; qh Þ; 0Þ; pffiffiffi
Proof. Following the proof of Theorem 1, we can show that
1 þ C 2pf
j2 1þ
C 2pf
ð22Þ
based on which we have the following similar error estimate:
r1
j2
ð21Þ
Proof. The proof is similar to that for Theorem 2 with some additional modifications. We omit the details. h
6 Ch
¼
sT k jDv h k20;T 2k jDv h k20;T
jku uh k1 þ jkp ph kdiv þ ska rðu uh Þk0
for all ðv h ; qh Þ 2 Vh Qh .
Fs ððv h ; qh Þ; 0Þ P
X
Applying (17)–(19), we obtain
ð18Þ
then we have the following coercivity estimate: there exists a constant C 3 > 0 independent of j and h such that
krv h k20 þ
jku uh k1 þ jkp ph kdiv þ hka rðu uh Þk0
8 T 2 Th :
Theorem 3. Assume that the assumption (18) holds. If for each element T in Th we take
sT ¼
1þ
C 2pf
1 ka rv h k20;T þ ka rv h k20;T 2 X j2 1 2 2 2 2 ¼ r v þ s j k D v k þ r v k ka k k T h 0 h 0;T h 0;T : 2 1 þ C 2pf T2Th
T2Th
Ls ððv h ; qh ÞÞ ¼ ðf ; r qh þ a rv h Þ0 þ
j2
krv h k20 þ
X
sT k jDv h þ a rv h k20;T
T2Th
krv h k20 þ
X
sT k jDv h k20;T
T2Th
þ 2ðjDv h ; a rv h Þ0;T þ ka rv h k20;T :
pffiffiffi pffiffiffi ððh þ j s þ h sÞkukrþ1 þ hkpkrþ1 Þ:
ð23Þ
We now give the numerical simulation for the simple one-dimensional problem (10) with the same setting as that used for the primitive LSFEM in Fig. 1. The results of the stabilized LSFEM with sT ¼ s :¼ 0:085 for all T 2 Th are depicted in Fig. 2. We observe that the boundary peaks seem to disappear, but the height of uh is lower than the primitive least-squares finite element solution. Motivated by this interesting observation, we perform the numerical simulation with sT ¼ s :¼ 0:085 for all T 2 Th . The approximate solution we obtain is presented in Fig. 3. Note that
187
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p stabilized LSFEM ph
1
1
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u stabilized LSFEM u
0.2
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 2. Exact solution vs. stabilized least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; s ¼ 0:085.
the peaks still appear near the boundaries, but away from the boundary the height of uh coincides with the exact solution very well. This observation is a motivation of the heuristic scheme called the negatively stabilized streamline diffusion LSFEM developed in Section 5.
Theorem 5. Assume that the assumption (18) holds. If for each element T in Th we take
dT ¼
e 2 Ch T
j
< 0;
then we have the following coercivity estimate: there exists a constant C 4 > 0 independent of j and h such that
4. The streamline diffusion LSFEM
j kv 2
2 h k1
j2 h2 2 þ kqh k2div þ h ka rv h k20 j2 þ j2 h2 þ h4
There is another way to smear out the peaks occurring near the boundaries for the primitive least-squares finite element solution. The underlying idea is similar to that for the streamline diffusion Galerkin method (see, e.g., [25]). We introduce another first-order system formulation for problem (1) which induces the so-called streamline diffusion (SD) least-squares finite element method [22,23]. Introducing the new dependent variable,
C4
p :¼ jru þ daðjDu þ a ruÞ;
Fd ððv h ; qh Þ; 0Þ
ð24Þ
Fd ððv h ; qh Þ; f Þ 2
¼ kr qh þ a rv h fd k0 þ
X
2
kqh þ jrv h dT aðjDv h þ a rv h Þk0;T ;
T2Th
Proof. Let a > 0 be a constant that will be determined later. For any ðv h ; qh Þ 2 Vh Qh , we have ¼ kr qh þ a rv h av h k20 þ 2aðr qh þ a rv h ; v h Þ0 a2 kv h k20 X þ kqh þ jrv h dT aðjDv h þ a rv h Þ arv h k20;T T2Th
X
þ
2aðqh þ jrv h dT aðjDv h þ a rv h Þ; rv h Þ0;T a2 krv h k20 :
T2Th
Employing the Green formula (3) and noting that ða rv h ; v h Þ0 ¼ 0, we obtain
Fd ððv h ; qh Þ; 0Þ P a2 kv h k20 þ 2ajkrv h k20
where the small parameter d which introduce a correction in the streamline direction is chosen element-wise and denoted by dT . Taking the variation, we obtain the so-called streamline diffusion LSFEM:
X
X
2adT ðjDv h ; a rv h Þ0;T
T2Th
2adT ka rv h k20;T a2 krv h k20 :
T2Th
Since for any positive constant c,
Determine ðuh ; ph Þ 2 Vh Qh such that Bd ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Ld ððv h ; qh ÞÞ 8 ðv h ; qh Þ 2 Vh Vh ;
ð28Þ
for all ðv h ; qh Þ 2 Vh Qh .
ð25Þ
We define the new least-squares energy functional Fd : Vh Qh ! R by
!
6 Fd ððv h ; qh Þ; 0Þ;
where d < 0 is an appropriate streamline diffusion parameter need to be determined later, then we have
r p þ a ru ¼ f þ dr ðf aÞ :¼ fd :
ð27Þ
ð26Þ
where the bilinear form Bd and linear form Ld are defined as follows, respectively:
Bd ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 X þ ðph þ jruh dT aðjDuh þ a ruh Þ; T2Th
X
P
Note that the streamline diffusion LSFEM is still a consistent scheme. For the streamline diffusion LSFEM (26), we have the following coercivity estimate:
X
1
adT j ck Dv h k20;T þ ka rv h k20;T ; c T2Th
we get that
Fd ððv h ; qh Þ; 0Þ
X
P 2ajkrv h k20 þ
qh þ jrv h dT aðjDv h þ a rv h ÞÞ0;T Ld ððv h ; qh ÞÞ ¼ ðfd ; r qh þ a rv h Þ0 :
2adT ðjDv h ; a rv h Þ0;T
T2Th
X T2Th
adT jckDv h k20;T
T2Th
adT 2
j ka rv h k20;T a2 ð1 þ C 2pf Þkrv h k20 ; c
where the Poincaré–Friedrichs inequality (4) has been used. Now j choosing a ¼ 2ð1þC 2 ; c ¼ j and using (27) with inverse estimate Þ pf (17), we obtain
188
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2 exact u negatively stabilized LSFEM u
0.2
exact p negatively stabilized LSFEM p
h
0
h
0
0
0.2
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 3. Exact solution vs. negatively stabilized least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; s ¼ 0:085.
Fd ððv h ; qh Þ; 0Þ P
j2 4ð1 þ C 2pf Þ
krv h k20 þ
X
e 2 Ch T
T2Th
2ð1 þ C 2pf Þ
ka rv h k20;T : ð29Þ
On the other hand,
kqh k20 6 2
X
X
kjrv h dT aðjDv h þ a rv h Þk20;T
T2Th
6 2Fd ððv h ; qh Þ; 0Þ þ 4j2 krv h k20 þ 8 þ8
X
j2 kv h k21 þ
jjdj jjdj þ 2d2 þ 1
! kqh k2div þ jjdjka rv h k20
6 Fs ððv h ; qh Þ; 0Þ;
ð33Þ
which implies the following error estimate:
kqh þ jrv h dT aðjDv h þ a rv h Þk20;T
T2Th
þ2
C4
X
d2T jaj2 j2 kDv h k20;T
T2Th
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi jjdj jku uh k1 þ kp ph kdiv þ jjdjka rðu uh Þk0 2 jjdj þ 2d þ 1 pffiffiffiffiffiffiffiffiffi r1 ð34Þ 6 Ch ðjjdj þ h jjdj þ hjdj þ hÞkukrþ1 þ hkpkrþ1 : In Fig. 4 we report the approximate solution ðuh ; ph Þ of the onedimensional problem (10) corresponding to the streamline diffusion LSFEM (23) with dT ¼ d :¼ 0:007 for all T 2 Th and the same values of j and h as before. Similar to the results obtained by stabilized LSFEM (16), boundary peaks disappear.
d2T jaj2 ka rv h k20;T :
T2Th
From (27) with inverse estimate (17), 2
kqh k20 6 2Fd ððv h ; qh Þ; 0Þ þ 4j2 krv h k20 þ Ch krv h k20 ! 2 2 h h þ C 2 Fd ððv h ; qh Þ; 0Þ ¼ C 1 þ 2 Fd ððv h ; qh Þ; 0Þ;
j
5. The negatively stabilized streamline diffusion LSFEM
ð30Þ
j
for some C > 0. Moreover, by the triangle inequality,
kr qh k20 6 2Fd ððv h ; qh Þ; 0Þ þ 2ka rv h k20 1 6 C 1 þ 2 Fd ððv h ; qh Þ; 0Þ: h Finally, the desired inequality (28) follows from (29)–(31).
ð31Þ
In the numerical results of streamline diffusion LSFEM reported in previous section, we find that the approximate solution is smooth but the height of uh is lower than the exact solution u. Motivated by the observation of the numerical results of the negatively stabilized LSFEM reported in Section 3, we now consider the streamline diffusion LSFEM with a negative stabilized term. Define the least-squares functional Fds : Vh Qh ! R by
Fds ððv h ; qh Þ; f Þ ¼ kr qh þ a rv h fd k20 þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2 h2 kp ph kdiv þ hka rðu uh Þk0 jku uh k1 þ j2 þ j2 h2 þ h4 ! 2 h r ð32Þ 6 Ch ð1 þ Þkukrþ1 þ kpkrþ1 :
j
Proof. The proof is similar to that for Theorem 2.h Remark 2. As in Remark 1, if Vh consists of piecewise linear polynomials ðr ¼ 1Þ and we take dT ¼ d for any d < 0 and for all T 2 Th , then the proof of Theorem 5 can be simplified to get the following estimate without assuming the inverse inequality (17):
kqh þ jrv h
T2Th
h
Theorem 6. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameter dT is defined as in (27), then the solution ðuh ; ph Þ of problem (26) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that
X
dT aðjDv h þ a rv h Þk20;T X þ sT k jDv h þ a rv h f k20;T ; T2Th
where dT < 0 and sT < 0 for all T 2 Th . Taking the variation, we get a negatively stabilized streamline diffusion LSFEM:
Determine ðuh ; ph Þ 2 Vh Qh such that Bds ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ Lds ððv h ; qh ÞÞ 8ðv h ; qh Þ 2 Vh Vh ;
ð35Þ
where the bilinear form Bds and linear form Lds are defined as follows, respectively:
Bds ððuh ; ph Þ; ðv h ; qh ÞÞ ¼ ðr ph þ a ruh ; r qh þ a rv h Þ0 X þ ðph þ jruh dT aðjDuh þ a ruh Þ; qh T2Th
þ jrv h dT aðjDv h þ a rv h ÞÞ0;T X þ sT ðjDuh þ a ruh ; jDv h þ a rv h Þ0;T ; T2Th
189
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact p SD LSFEM p
1
1
h
0.8
0.8
0.6 0.6 0.4 0.4
0.2 exact u SD LSFEM uh
0.2 0
0
0.2
0
0.4
0.6
0.8
1
−0.2
0
0.2
0.4
0.6
x
0.8
1
x
Fig. 4. Exact solution vs. streamline diffusion least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; d ¼ 0:007.
Lds ððv h ; qh ÞÞ ¼ ðfd ; r qh þ a rv h Þ0 þ
X
sT ðf ; jDv h þ a rv h Þ0;T :
Proof. The proof is similar to that for Theorem 6. h
T2Th
Still, the negatively stabilized streamline diffusion LSFEM (35) is a consistent scheme and we have the following estimate: Theorem 7. Assume that the assumption (17) holds and for each element T in Th the parameter dT is defined as in (27). If take
sT ¼
e 2 Ch T eÞ 8ð1 þ C 2pf Þð1 þ C
< 0;
ð36Þ
then we have the following coercivity estimate: there exists a constant C 5 > 0 independent of j and h such that
C5
j2 kv h k21 þ
j2 h2 2 kqh k2div þ h ka rv h k20 2 j þ j2 h2 þ h4
!
6 Fds ððv h ; qh Þ; 0Þ;
ð37Þ
C5
j2 kv h k21 þ
!
jjdj jjdj þ 2d2 þ 1
kqh k2div þ jjdjka rv h k20
6 Fds ððv h ; qh Þ; 0Þ:
ð40Þ
Based on (40), we can derive the error estimate: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi jjdj kp ph kdiv þ jjdjka rðu uh Þk0 jku uh k1 þ 2 jjdj þ 2d þ 1 qffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi r1 6 Ch ð jjdj þ jdj þ 1Þh þ jjdj þ j3 jdj kukrþ1 þ hkpkrþ1 :
ð41Þ
for all ðv h ; qh Þ 2 Vh Qh . Proof. Similar to the proofs of Theorems 3 and 5, one can derive that
Fds ððv h ; qh Þ; 0Þ P
Remark 3. Assume that Vh consists of piecewise linear polynomials ðr ¼ 1Þ. If for each element T in Th we take dT ¼ d for any d < 0 jd < 0, then we have the following estimate withand sT ¼ s :¼ 1þC 2 pf out assuming the inverse inequality (17):
j2 C 2pf Þ
4ð1 þ X þ2
krv h k20 þ
X
e h2 C T
T2Th
2ð1 þ C 2pf Þ
ka rv h k20;T
sT ðk jDv h k20;T þ ka rv h k20;T Þ:
In Fig. 5 it is shown the negatively stabilized streamline diffusion least-squares finite element solution with dT ¼ d :¼ 0:007 and sT ¼ s :¼ 0:085 for all T 2 Th . Although the magnitude of s is larger than our estimate given in Remark 3, but it seems to achieve the optimal result for mesh size h ¼ 1=32. 6. Numerical experiments for 2D problems
T2Th
From (36), we conclude that
Fds ððv h ; qh Þ; 0Þ P
j2 8ð1 þ
C 2pf Þ
krv h k20 þ
X
e h2 C T
T2Th
4ð1 þ C 2pf Þ
ka rv h k20;T :
In this section, we perform a series of numerical simulations for four test problems by using the developed methods, including the primitive LSFEM (6), stabilized LSFEM (16), streamline diffusion LSFEM (26), and the negatively stabilized streamline diffusion LSFEM (35).
ð38Þ Moreover, we still have (30) and (31) with Fd ððv h ; qh Þ; 0Þ replaced by Fds ððv h ; qh Þ; 0Þ, and then the desired estimate (37) follows from (38), (30) and (31). h Similarly, we have the following error estimate: Theorem 8. Assume that problem (2) has a sufficiently smooth solution ðu; pÞ 2 ðH10 ðXÞ Hðdiv; XÞÞ \ ðHrþ1 ðXÞÞ3 . If the assumption (18) holds and for each element T in Th the parameters dT is defined as in (27) and sT as in (36), then the solution ðuh ; ph Þ of problem (35) satisfies the following error estimate: there exists a constant C > 0 independent of j and h such that
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2 h2 jku uh k1 þ kp ph kdiv þ hka rðu uh Þk0 j2!þ j2 h2 þ h4 ! 2 h r kukrþ1 þ kpkrþ1 : ð39Þ 6 Ch 1þ
j
6.1. An example with analytic solution We first consider a simple problem with a ¼ ð0; 1Þ> on the unit square X ¼ ð0; 1Þ ð0; 1Þ subjected to the following boundary conditions [8]:
8 > < x ¼ 0 and 0 6 y 6 1; u ¼ 0 when x ¼ 1 and 0 6 y 6 1; > : y ¼ 1 and 0 6 x 6 1;
u ¼ sinðpxÞ when y ¼ 0 and 0 6 x 6 1: One can check that the exact solution is given by
u¼
1 ðem2 m1 em1 y em2 y Þ sinðpxÞ; em2 m1 1
where
m1 ¼
1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4j2 p2 2j
and m2 ¼
1þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4j2 p2 : 2j
190
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2 exact u negatively stabilized SD LSFEM uh
0.2 0
0
0.2
0.4
0.6
0.8
exact p negatively stabilized SD LSFEM p
h
0 1
−0.2
0
0.2
x
0.4
0.6
0.8
1
x
Fig. 5. Exact solution vs. negatively stabilized streamline diffusion least-squares finite element solution: h ¼ 1=32; j ¼ 102 ; d ¼ 0:007; s ¼ 0:085.
Table 1 Relative error behavior for
j ¼ 1 and r ¼ 1.
Method
u=p
h ¼ 1=32
h ¼ 1=64
h ¼ 1=128
h ¼ 1=256
LSFEM
u
2
L
9.3538E4
2.3849E4
6.0488E5
1.5344E5
1.98
u
H1
2.0763E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9168E2
7.2544E3
2.7362E3
1.0290E3
1.41
L2
4.9617E2 9.1326E4
2.5847E2 2.3295E4
1.3219E2 5.9100E5
6.6925E3 1.5001E5
0.96 1.98
u
H1
2.0762E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9167E2
7.2544E3
2.7362E3
1.0290E3
1.41
L2
4.9617E2 9.2450E4
2.5847E2 2.3577E4
1.3219E2 5.9807E5
6.6925E3 1.5176E5
0.96 1.98
u
H1
2.0762E2
1.0379E2
5.1894E3
2.5947E3
1.00
p
L2 HðdivÞ
1.9166E2
7.2542E3
2.7361E3
1.0290E3
1.41
L2
4.9617E2 9.3343E4
2.5847E2 2.3801E4
1.3219E2 6.0372E5
6.6925E3 1.5314E5
0.96 1.98
p u
Stabilized LSFEM
p u
SD LSFEM
p u
Negatively stabilized
Norm
Order
u
H1
2.0763E2
1.0379E2
5.1895E3
2.5947E3
1.00
SD
p
1.9166E2
7.2541E3
2.7360E3
1.0290E3
1.41
LSFEM
p
L2 HðdivÞ
4.9617E2
2.5847E2
1.3219E2
6.6925E3
0.96
(Stabilized LSFEM:
2
2
2
2
s ¼ 0:1h ; SD LSFEM: d ¼ 0:05h =j; negatively stabilized SD LSFEM: d ¼ 0:07h =j; s ¼ 0:06h ).
Table 2 Relative error behavior for
j ¼ 102 and r ¼ 1.
Method
u=p
Norm
h ¼ 1=32
h ¼ 1=64
h ¼ 1=128
h ¼ 1=256
Order
LSFEM
u
L2
9.7176E2
3.0693E2
8.2541E3
2.1150E3
1.84
u
H1
4.4607E5
2.9571E5
1.6109E5
8.2399E2
0.81
p
L2 HðdivÞ
3.8496E5
1.2947E5
3.6118E2
9.7787E3
1.77
L2
4.9399E5 1.7826E5
3.2868E5 7.0424E2
1.7922E5 2.0870E2
9.1710E2 5.4864E3
0.81 1.67
u
H1
4.3500E5
2.9518E5
1.6126E5
8.2437E2
0.80
p
L2 HðdivÞ
3.9890E5
1.3585E5
3.8033E2
1.0299E2
1.76
L2
4.7611E5 6.8719E2
3.2665E5 2.7989E2
1.7917E5 8.3983E3
9.1720E2 2.2232E3
0.79 1.65
u
H1
4.1448E5
2.8767E5
1.5982E5
8.2231E2
0.78
p
L2 HðdivÞ
5.1168E5
1.6527E5
4.5548E2
1.2106E2
1.80
L2
4.5991E5 4.3060E2
3.1980E5 1.3619E2
1.7781E5 3.6980E3
9.1522E2 9.6171E4
0.78 1.83
u
H1
4.1941E5
2.8944E5
1.6012E5
8.2270E2
0.78
p
L2 HðdivÞ
5.8724E5
1.8574E5
5.0631E2
1.3325E2
1.82
4.6612E5
3.2198E5
1.7818E5
9.1571E2
0.78
p u
Stabilized LSFEM
p u
SD LSFEM
Negatively stabilized SD LSFEM (Stabilized LSFEM:
p u
p
s ¼ 0:1h2 ; SD LSFEM: d ¼ 0:05h2 =j; negatively stabilized SD LSFEM: d ¼ 0:07h2 =j; s ¼ 0:06h2 ).
Substituting the above exact solution into problem (2), we can readily get the right-hand side function f. We consider two different values for the diffusivity, namely, j ¼ 1 and j ¼ 102 . For simplifying the numerical implementation, we shall assume that the unit square domain X is uniformly 2 partitioned into a set of h square subdomains with side-length
h ¼ 2‘ for ‘ ¼ 5; 6; 7; 8. We will apply continuous piecewise bilinear finite elements ðr ¼ 1Þ to approximate all components of the exact solution for all least-squares finite element schemes. In all computations, a double precision conjugate gradient solver is applied to solve the associated linear systems. We compute the relative error in the L2 ðXÞ-norm and H1 ðXÞ-norm for u, and L2 ðXÞ-norm
191
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
exact solution
exact solution 1 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0
y
0.5
x
0
1 1
0
0.2
Fig. 6. Exact solution u and its contours:
LSFEM
0.4
0.6
0.8
1
j ¼ 102 .
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
1 1
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 7. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y 1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 8. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; s ¼ 0:1h .
SD LSFEM
SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 9. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; d ¼ 0:05h =j.
and Hðdiv; XÞ-norm for p. The results are collected in Tables 1 and 2. The elevation and contour plots for the exact solution u and
approximate solutions uh for h ¼ 1=32 and in Figs. 6–10.
j ¼ 102 are displayed
192
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y 1 1
0.5
x
0
0
0.2
0.4
0.6
0.8
1 2
2
Fig. 10. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 102 ; d ¼ 0:07h =j; s ¼ 0:06h .
LSFEM
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0
0.5 y 1 1
x
0.5
0
0
0.2
0.4
0.6
0.8
1 6
Fig. 11. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x 0.5
0
0
0.2
0.4
0.6
0.8
1 2
Fig. 12. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 106 ; s ¼ 0:2h .
SD LSFEM
SD LSFEM
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
1 1
x 0.5
0
0
0.2
0.4
0.6
0.8
1 6
2
Fig. 13. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 ; d ¼ 0:00002h =j.
Numerical results reported in Table 1 for j ¼ 1 confirm the theoretical error estimates with r ¼ 1. That is, all the primitive LSFEM and its variants developed in this paper achieve optimal conver-
gence in the H1 ðXÞ-norm for u and in the Hðdiv; XÞ-norm for p. Furthermore, an examination of the numerical results given in Table 1 shows that the convergence rate for u in L2 ðXÞ-norm seems to be
193
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
1 1
0.5
0
0
0.2
0.4
0.6
0.8
1 2
2
Fig. 14. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 106 ; d ¼ 0:0000163h =j; s ¼ 0:06h .
LSFEM
LSFEM 1 0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
1 0.8 0.6 0.4 0.2 0
0 0 0.5
y
0.5
x
0
1 1
0
0.2
0.4
0.6
0.8
1 4
Fig. 15. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 10 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
1 0.8 0.6 0.4 0.2 0
0 0 0.5
y
0.5
x
0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 16. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; s ¼ 1.
SD LSFEM
SD LSFEM 1
1 0.8
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.6 0.4 0.2 0
0 0 0.5
0.5 y
1 1
x
0
0 0
0.2
0.4
0.6
0.8
1
Fig. 17. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 104 .
optimal, although we have not provided the analysis. Also, from Table 2 for j ¼ 102 , it is no surprise that the convergence behavior
for all unknowns is getting worse when the diffusivity j is getting small.
194
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1
1 0.8
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.6 0.4 0.2 0
0 0 0.5
y
x
0.5 0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 18. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 104 ; s ¼ 0:5.
LSFEM
LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
0.5
0
1 1
0 0
0.2
0.4
0.6
0.8
1
Fig. 19. Least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 .
stabilized LSFEM
stabilized LSFEM 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.8 0.6 0.4 0.2 0 0 0.5
0 y
x
0.5
1 1
0
0 0
0.2
0.4
0.6
0.8
1 2
Fig. 20. Stabilized least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; s ¼ 0:2h .
6.2. Boundary layer problem
6.3. Interior layer problem
We consider the test case where X ¼ ð0; 1Þ ð0; 1Þ; j ¼ 106 ; a ¼ ð0; 1Þ> , f ¼ 0, and the boundary conditions are given by
ary data with a uniform velocity field a ¼
u ¼ 1 on fðx; yÞ : 0 6 x 6 1; y ¼ 0g; u ¼ 0 on fðx; yÞ : 0 6 x 6 1; y ¼ 1g; and the homogeneous Neumann boundary condition are imposed on the other two sides, namely,
p1 ¼ 0 on fðx; yÞ : x ¼ 0; 1 and 0 6 y 6 1g: A uniform mesh with h ¼ 25 and piecewise bilinear finite elements are employed for all least-squares methods. The results are shown in Figs. 11–14. One can observe that the negatively stabilized streamline diffusion least-squares method is able to better capture the boundary layer behavior when compared with other leastsquares methods.
Let us consider the propagation of a discontinuity in the bound > p1ffiffi ; p1ffiffi on the unit 2 2
square domain X ¼ ð0; 1Þ ð0; 1Þ with f ¼ 0 and boundary conditions are given by
j ¼ 104 . The
u ¼ 1 on fðx; yÞ : x ¼ 0; 0 6 y 6 1g [ fðx; yÞ : 0 6 x 6 1; y ¼ 1g; u ¼ 0 on fðx; yÞ : 0 < x 6 1; y ¼ 0g [ fðx; yÞ : x ¼ 1; 0 6 y < 1g: We have an interior layer due to the discontinuity of the boundary data. A uniform mesh with h ¼ 25 and piecewise bilinear finite elements are employed for all least-squares methods. In Figs. 15–18, the elevation and contour plots for approximate solution uh produced by various methods are displayed. Except the primitive LSFEM, all other least-squares methods perform very well.
195
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
SD LSFEM
SD LSFEM 1
1
0.8
0.8
0.6
0.6
0.8 0.6 0.4
0.4
0.4
0.2 0 0
0.2
0.2 0.5
0 y
x
0.5
0
1 1
0 0
0.2
0.4
0.6
0.8
1 2
Fig. 21. Streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 0:0017h =j.
negatively stabilized SD LSFEM
negatively stabilized SD LSFEM 1 1
1
0.8 0.8
0.8 0.6 0.4
0.6
0.6
0.4
0.4
0.2
0.2
0.2 0 0 0.5
0 y 1 1
x
0.5
0
0 0
0.2
0.4
0.6
0.8
1 2
2
Fig. 22. Negatively stabilized streamline diffusion least-squares finite element solution uh and its contours: h ¼ 1=32; j ¼ 104 ; d ¼ 0:0017h =j; s ¼ 0:055h .
6.4. Boundary and interior layer problem We now construct a model problem possessing both boundary and interior layer behavior on the unit square domain X ¼ ð0; 1Þ ð0; 1Þ. As in Section 6.3, we consider the propagation of a discontinuity in the boundary data with constant velocity field pffiffi> a ¼ 12 ; 23 of size one forming a 60° with the x-axis. In this example, f ¼ 0; j ¼ 104 , and the boundary conditions are given by
the primitive LSFEM performs poorly for the convection-dominated model while the other three variants perform considerably better for interior layer problems, and the negatively stabilized streamline diffusion LSFEM is able to better capture the boundary layer behavior when compared with other LSFEMs. But for problem possessing both interior and boundary layer structures in the solution, more investigation is needed in order to reach an efficient LSFEM. The study in this direction is in progress.
u ¼ 1 on fðx; yÞ : x ¼ 0; 0 6 y < 1g [ fðx; yÞ : 0 6 x 6 1=2; y ¼ 0g; u ¼ 0 on fðx; yÞ : 0 6 x 6 1; y ¼ 1g [ fðx; yÞ : x ¼ 1; 0 6 y 6 1g
Acknowledgment
[ fðx; yÞ : 1=2 < x 6 1; y ¼ 0g:
This work was partially supported by the National Science Council of Taiwan.
The elevation and contour plots of all approximate solutions are displayed in Figs. 19–22. Unfortunately, all the primitive least-squares method and its variants do not give reasonable results for problem possessing both interior and boundary layer structures in the solution. 7. Concluding remarks In this paper, we have proposed and analyzed the primitive LSFEM and its three variants, including the stabilized LSFEM, streamline diffusion LSFEM and the negatively stabilized streamline diffusion LSFEM. The coercivity (stability) estimates for all methods are given. One of the common advantages of these least-squares finite element methods is that the resulting linear system is symmetric and positive definite. It was observed that
References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convective dominated flows with a particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Engrg. 32 (1982) 199–259. [2] L.P. Franca, S.L. Frey, T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective–diffusive model, Comput. Meth. Appl. Mech. Engrg. 95 (1992) 253–276. [3] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations, Comput. Meth. Appl. Mech. Engrg. 73 (1989) 173–189. [4] L.P. Franca, G. Hauke, A. Masud, Revisiting stabilized finite element methods for the advective–diffusive equation, Comput. Meth. Appl. Mech. Engrg. 195 (2006) 1560–1572. [5] F. Brezzi, A. Russo, Choosing bubbles for advection–diffusion problems, Math. Mod. Meth. Appl. Sci. 4 (1994) 571–587.
196
P.-W. Hsieh, S.-Y. Yang / Comput. Methods Appl. Mech. Engrg. 199 (2009) 183–196
[6] F. Brezzi, L.P. Franca, A. Russo, Further considerations on residual-free bubbles for advective–diffusive equations, Comput. Meth. Appl. Mech. Engrg. 166 (1998) 25–33. [7] F. Brezzi, T.J.R. Hughes, L.D. Marini, A. Russo, E. Süli, A priori error analysis of residual-free bubbles for advection–diffusion problems, SIAM J. Numer. Anal. 36 (1999) 1933–1948. [8] L.P. Franca, A. Nesliturk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by a two-level finite element method, Comput. Meth. Appl. Mech. Engrg. 166 (1998) 35–49. [9] P.B. Bochev, M.D. Gunzburger, Finite element methods of least-squares type, SIAM Rev. 40 (1998) 789–837. [10] B.-N. Jiang, The Least-Squares Finite Element Method, Springer-Verlag, Berlin, 1998. [11] P.B. Bochev, M.D. Gunzburger, Analysis of least-squares finite element methods for the Stokes equations, Math. Comp. 63 (1994) 479–506. [12] Z. Cai, R. Lazarov, T.A. Manteuffel, S.F. McCormick, First-order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal. 31 (1994) 1785–1799. [13] Z. Cai, T.A. Manteuffel, S.F. McCormick, First-order system least squares for second-order partial differential equations: Part II, SIAM J. Numer. Anal. 34 (1997) 425–454. [14] C.L. Chang, B.-N. Jiang, An error analysis of least-squares finite element method of velocity–pressure–vorticity formulation for Stokes problem, Comput. Meth. Appl. Mech. Engrg. 84 (1990) 247–255. [15] C.L. Chang, S.-Y. Yang, C.-H. Hsu, A least-squares finite element method for incompressible flow in stress–velocity–pressure version, Comput. Meth. Appl. Mech. Engrg. 128 (1995) 1–9.
[16] C.L. Chang, S.-Y. Yang, Analysis of the L2 least-squares finite element method for the velocity–vorticity–pressure Stokes equations with velocity boundary conditions, Appl. Math. Comput. 130 (2002) 121–144. [17] J.M. Deang, M.D. Gunzburger, Issues related to least-squares finite element methods for the Stokes equations, SIAM J. Sci. Comput. 20 (1998) 878–906. [18] A.I. Pehlivanov, G.F. Carey, R.D. Lazarov, Least-squares mixed finite elements for second-order elliptic problems, SIAM J. Numer. Anal. 31 (1994) 1368–1377. [19] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994. [20] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, SpringerVerlag, New York, 1991. [21] J.M. Fiard, T.A. Manteuffel, S.F. McCormick, First-order system least squares (FOSLS) for convection–diffusion problems: numerical results, SIAM J. Sci. Comput. 19 (1998) 1958–1979. [22] R.D. Lazarov, L. Tobiska, P.S. Vassilevski, Streamline diffusion least-squares mixed finite element methods for convection–diffusion problems, East-West J. Numer. Math. 5 (1997) 249–264. [23] R.D. Lazarov, P.S. Vassilevski, Least-squares streamline diffusion finite element approximations to singularly perturbed convection–diffusion problems, in: L.G. Vulkov, J.J.H. Miller, G.I. Shishkin (Eds.), Analytical and Numerical Methods for Singularly Perturbed Problems, Nova Science Publishing House, 2000, pp. 83–94. [24] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, New York, 1986. [25] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, Cambridge, 1987.