- Email: [email protected]
Mixed finite element method and high-order local artificial boundary conditions for exterior problems of elliptic equation
Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027 www.elsevier.com/locate/cma
Mixed finite element method and high-order local artificial boundary conditions for exterior problems of elliptic equation q Houde Han *, Chunxiong Zheng Department of Mathematical Sciences, Tsinghua University, Beijing 100084, People’s Republic of China Received 21 February 2001; received in revised form 18 October 2001; accepted 19 October 2001
Abstract The mixed finite element method is used to solve the exterior elliptic problem with high-order local artificial boundary conditions. New unknowns are introduced to reduce the order of the derivatives to two. This leads to an equivalent mixed variational problem such that the normal finite element can be used and special finite elements are no longer needed on the adjacent layer of the artificial boundary. Error estimates are obtained for some local artificial boundary conditions with prescribed order. Numerical examples are presented and the results demonstrate the effectiveness of this method. 2001 Elsevier Science B.V. All rights reserved. Keywords: Mixed finite element method; Artificial boundary condition; Elliptic equation
1. Introduction Elliptic problems on exterior domain arise in many application fields. For example, the incompressible irrotational flow around a body is described by the exterior problems of Laplace equation. There are many numerical methods to solve these exterior problems. The most popular is the artificial boundary method. In the procedure of this method, an artificial boundary is firstly introduced to divide the exterior domain into a bounded part and an unbounded part. After a suitable artificial boundary condition is set up on the artificial boundary, the original exterior problem is reduced to a new one only defined on the bounded domain, and then it can be solved by a proper numerical method. In the last two decades, many authors have worked on this subject for various problems by different techniques (see [7–9,12–25] and the references
q
This work was supported partly by the Special Funds for Major State Basic Research Projects of China and the National Natural Science Foundation of China. Computation was supported by the State Key Laboratory of the Scientific and Engineering Computing in China. * Corresponding author. E-mail address: [email protected] (H. Han). 0045-7825/01/$ - see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 3 6 5 - 6
2012
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
therein). Generally speaking, the artificial boundary conditions can be classified into nonlocal artificial boundary conditions, local artificial boundary conditions and discrete artificial boundary conditions. A thorough review concerning nonlocal and local artificial boundary conditions can be found in [12]. For the discrete artificial boundary conditions, the readers can refer to [2,3,13,18–22]. For high-order local artificial boundary conditions, since the appearance of high-order derivatives, the new difficulty is posed in the implementation with a conforming finite element method. To allow the use of high-order local artificial boundary conditions, Givoli and others [10,12] developed a family of special elements and in [11], Givoli and Patlashenko proposed a new finite element scheme for the numerical simulation of time-harmonic wave scattering problems in unbounded domains. In this paper, mixed finite element method will be used to overcome this new difficulty.
2. High-order local artificial boundary conditions Suppose Ci is a bounded, simple closed curve in R2 and X is the exterior domain with boundary Ci . The boundary Ci is decomposed into two disjoint parts: CD with a nonzero measure where a Dirichlet boundary condition is given and CN where a Neumann boundary condition is given. Consider the following exterior problems of a linear elliptic second-order equation: r ðjðxÞruðxÞÞ þ bðxÞuðxÞ ¼ f ðxÞ in X;
ð1Þ
u¼g
ð2Þ
on CD ;
ou ¼k on
ð3Þ
on CN ;
u is bounded
when jxj ! 1;
ð4Þ
, g and k are the where uðxÞ is the unknown function, jðxÞ P 1, bðxÞ P 0, f ðxÞ are the given functions on X given functions on CD and CN . Furthermore, suppose suppðf Þ, suppðj 1Þ and suppðbÞ are compact, jðxÞ 2 L1 ðXÞ, bðxÞ 2 L1 ðXÞ, f ðxÞ 2 L2 ðXÞ, g 2 H 3=2 ðCD Þ, kðxÞ 2 H 1=2 ðCN Þ. For brevity let gðxÞ 0 in the following. Let R0 ¼ 0 þ maxfr: ðr; hÞ 2 Ci [ suppðf Þ [ suppðjÞ [ suppðbÞg; where 0 is a positive constant and C0 ¼ fðR0 ; hÞ: 0 6 h 6 2pg. Then C0 X. Introduce a circular artificial boundary Ce fðR; hÞ: 0 6 h 6 2pg with R > R0 . Ce divides the exterior domain into two parts: the e . On the unbounded unbounded part Xe ¼ fðr; hÞ: r > R; 0 6 h 6 2pg and the bounded part Xi ¼ X n X domain Xe , jðxÞ 1, bðxÞ 0, f ðxÞ 0 and Eq. (1) is reduced to a Laplace equation. Now recall the highorder local artificial boundary conditions on the given artificial boundary Ce [4,12,23]. Suppose uðxÞ is the solution of problem (1)–(4). Then on X0 fðr; hÞ: r > R0 ; 0 6 h 6 2pg, uðxÞ satisfies Laplace equation and condition (4). By separation of variables, uðxÞ can be given on X0 by n 1 a0 X R0 ðan cos nh þ bn sin nhÞ ð5Þ uðr; hÞ ¼ þ 2 r n¼1 with 1 an ¼ p
Z 0
2p
1 uðR0 ; hÞ cos nh dh ¼ p
R R0
n Z
2p
uðR; hÞ cos nh dh; 0
ð6Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
bn ¼
1 p
Z
2p
uðR0 ; hÞ sin nh dh ¼ 0
1 p
R R0
n Z
2013
2p
ð7Þ
uðR; hÞ sin nh dh: 0
Therefore n 1 ou ouðR; hÞ 1X R0 ¼ n ðan cos nh þ bn sin nhÞ on Ce or R n¼1 R Z 2p 1 1 X n uðR; h0 Þ cos nðh h0 Þ dh0 : ¼ pR n¼1 0
ð8Þ
This is the desired exact boundary condition on Ce . Furthermore, n 1 a0 X R0 ðan cos nh þ bn sin nhÞ uðR; hÞ ¼ þ 2 R n¼1
ð9Þ
and n 1 o2m uðR; hÞ X 2 m R0 ¼ ðn Þ ðan cos nh þ bn sin nhÞ R oh2m n¼1
ð10Þ
with positive integer m. A computation shows
N N 2m ouðR; hÞ 1 X o2m uðR; hÞ ouðR; hÞ 1 X m m N o uðR; hÞ ð1Þ aNm ¼ ð1Þ a m on R m¼1 or R m¼1 oh2m oh2m n 1 1X R0 ¼ An ðan cos nh þ bn sin nhÞ; R n¼1 R
where N is a positive integer and faN1 ; aN2 ; . . . ; aNN g are N constants and An ¼ n
N X
aNm n2m ;
n ¼ 1; 2; . . .
m¼1
Let faN1 ; aN2 ; . . . ; aNN g be the solution of the following linear algebraic equations: N X
aNm n2m ¼ n;
n ¼ 1; 2; . . . ; N :
ð11Þ
m¼1
Then
n N 1 ouðR; hÞ 1 X o2m uðR; hÞ 1 X R0 1 m ð1Þ aNm ¼ A ða cos nh þ b sin nhÞ ¼ O : n n n on R m¼1 R n¼1 RN þ2 R oh2m
Table 1 Solution of algebraic equations (11) aN1 N ¼1 N ¼2 N ¼3 N ¼4 N ¼5
aN2
aN3
aN4
1 60 11 360 71 1728
1 1008
aN5
1 7 6 37 30 533 420 1627 1260
16 14 43 144 107 324
13 6048
1 25 920
ð12Þ
2014
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Omitting the RHS of (12) gives the following high-order local artificial boundary conditions on the given artificial boundary, which are exact for the first N modes: N 2m ou 1X m N o uðR; hÞ ¼ ð1Þ a : ð13Þ m on Ce R m¼1 oh2m For any positive integer N, it is easy to be verified that algebraic equations (11) have a unique solution. The solution for N ¼ 1; 2; . . . ; 5 is illustrated in Table 1. Applying the local artificial boundary conditions gives a series of approximate boundary value problems defined on the bounded domain Xi : r ðjðxÞruN ðxÞÞ þ bðxÞuN ðxÞ ¼ f ðxÞ
ð14Þ
in Xi ;
uN ¼ 0 on CD ;
ð15Þ
ouN ¼k on
ð16Þ
on CN ;
N ouN 1X o2m uN m ¼ ð1Þ aNm R m¼1 on oh2m
ð17Þ
on Ce :
Let H m ðXi Þ and H s ðCe Þ be the usual Sobolev spaces on Xi and Ce with integer m and real number s [1]. Suppose H1 ðXi Þ ¼ fv 2 H 1 ðXi Þ: v ¼ 0 on CD g: Then the weak form of problem (1)–(3) with the exact boundary condition (8) is: Find u 2 H1 ðXi Þ such that aðu; vÞ þ bðu; vÞ ¼ f ðvÞ 8v 2 H1 ðXi Þ; where aðu; vÞ ¼
Z
jðxÞruðxÞ rvðxÞ dx þ
Xi
Z
ð18Þ
bðxÞuðxÞvðxÞ dx 8u; v 2 H 1 ðXi Þ;
Z Z 1 X n 2p 2p cos nðh h0 ÞuðR; hÞvðR; h0 Þ dh dh0 bðu; vÞ ¼ p 0 0 n¼1 f ðvÞ ¼
Z
ð19Þ
Xi
f ðxÞvðxÞ dx þ
Xi
Z
8u; v 2 H 1 ðXi Þ;
jðxÞkðxÞvðxÞ ds 8v 2 H 1 ðXi Þ:
ð20Þ
ð21Þ
CN
In order to formulate the weak form of the problem (14)–(17), define VN ¼ fv 2 H1 ðXi Þ: vjCe 2 H N ðCe Þg with norm kvkVN jvj1;Xi þ jvjN ;Ce :
ð22Þ
Let bN ðu; vÞ ¼
N X m¼1
aNm
Z
2p 0
om uðR; hÞ om vðR; hÞ dh: ohm ohm
ð23Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2015
Then the boundary value problem (14)–(17) is equivalent to the following variational problem: Find uN 2 VN such that aðuN ; vÞ þ bN ðuN ; vÞ ¼ f ðvÞ
8v 2 VN :
ð24Þ
Theorem 1. Suppose N is an odd number: • The bilinear form að; Þ þ bN ð; Þ is coercive on VN and the variational problem (24) is well-posed. • Let uðxÞ be the solution of problem (1)–(4). Then the restriction of uðxÞ on Xi is the solution of variational problem (18). Assume that uN is the solution of the problem (24). Then the following error estimate holds: ku uN kVN 6 CN
R0 R
N þ1 jujN ;C0 ;
ð25Þ
where CN is a constant only depending on N.
Remark. In [4], this theorem is proved under the hypothesis of 1 6 N < 20. But the restriction on N is not essential. The key point is to prove the coercivity property of bN ð; Þ on H N ðCe Þ, which has been done in [26] for all odd numbers. So the proof of this theorem is omitted here. When N > 1, the high-order derivatives of the unknown function on boundary Ce appear in the variational problem (24), therefore, the standard C ð0Þ -element cannot be used to get the numerical approximation of problem (24). Givoli and others [10,12] developed a family of special elements and their extension, a hierarchy of finite elements, to overcome this new difficulty. In the following sections, a different approach, mixed finite element method, will be proposed to compute the numerical solution of the boundary value problem (14)–(17) for N > 1.
3. Mixed variational formulation of problem (14)–(17) Suppose that N ¼ 2L þ 1 is an odd number and J ðwÞ ¼ 12ðaðw; wÞ þ bN ðw; wÞÞ f ðwÞ
8w 2 VN :
According to Theorem 1, the boundary value problem (14)–(17) is equivalent to the following minimization problem: Find uN 2 VN such that J ðuN Þ ¼ inf J ðwÞ:
ð26Þ
w2VN
For 8w 2 VN , introduce new unknowns on the artificial boundary Ce denoted by wk ðR; hÞ ¼
o2k wðR; hÞ ; oh2k
k ¼ 1; 2; . . . ; L
ð27Þ
and w0 ðr; hÞ ¼ wðr; hÞ
8ðr; hÞ 2 Xi :
ð28Þ
2016
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
The bilinear form bN ðw; wÞ can be rewritten as follows: 2 Z 2p m N X o wðR; hÞ N bN ðw; wÞ ¼ am dh ohm 0 m¼1 2 2 Z 2p 2kþ1 Z 2p 2k L L X X o wðR; hÞ o wðR; hÞ N N ¼ a2kþ1 dh þ a2k dh oh2kþ1 oh2k 0 0 k¼0 k¼1 2 Z 2p Z 2p L L X X owk ðR; hÞ N ðwk ðR; hÞÞ2 dh: ¼ a2kþ1 dh þ aN2k oh 0 0 k¼0 k¼1 Denote H01 ðCe Þ ¼
w 2 H 1 ðCe Þ and
Z
w¼0 ;
Ce
n o ¼ ðw0 ; w1 ; w2 ; . . . ; wL Þ j w0 2 V1 ; wk 2 H01 ðCe Þ; k ¼ 1; 2; . . . ; L ; WN ¼ w n o ¼ ðq1 ; q2 ; . . . ; qL Þ j qk 2 H01 ðCe Þ; k ¼ 1; 2; . . . ; L Q¼ q and Þ ¼ bN ð u; w
L X
aN2kþ1
Z 0
k¼0
2p
L X ouk ðR; hÞ owk ðR; hÞ dh þ aN2k oh oh k¼1
Z
2p
uk ðR; hÞwk ðR; hÞ dh;
0
Þ ¼ aðu0 ; w0 Þ þ bN ð Þ; a ð u; w u; w f ð wÞ ¼ f ðw0 Þ; Þ f ð J ð wÞ ¼ 12a ð w; w wÞ; Cð w; qÞ ¼
L Z X k¼1
2p 0
owk1 oqk w k qk þ dh oh oh
8 w 2 WN ; 8q 2 Q;
2 WN and q 2 Q. Then where u 2 WN , w Þ w; w bN ðw; wÞ ¼ bN ð and Eq. (27) are equivalent to the following variational form Cð w; qÞ ¼ 0
8 q 2 Q:
ð29Þ
Furthermore introduce n o jw 2 WN and Cð w; qÞ ¼ 0; 8 q2Q : SN ¼ w Obviously SN is a subspace of WN and the minimization problem (26) is equivalent to the following conditional minimization problem: Find u 2 SN such that uÞ ¼ min J ð wÞ: J ð 2SN w
ð30Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2017
Introduce the Lagrange multiplier q 2 Q and let Þ þ Cð Ið w; qÞ ¼ 12a ð w; w w; qÞ f ð wÞ:
ð31Þ
Then the conditional minimization problem (30) is equivalent to the following saddle problem [5]: Find ð uN ; pN Þ 2 WN Q such that pN Þ ¼ inf sup Ið w; qÞ: Ið uN ;
ð32Þ
2WN w q2Q
Furthermore the saddle problem is equivalent to the following mixed variational problem: Find ð uN ; pN Þ 2 WN Q such that Þ þ Cð a ð uN ; w w; pN Þ ¼ f ð wÞ 8 w 2 WN ; Þ ¼ 0 8 Cð uN ; q q 2 Q:
ð33Þ
Finally the mixed variational problem (33) is obtained, which is equivalent to the boundary value problem (14)–(17). And now there are no high-order derivatives of the unknowns involved in the formulation. Lemma 1. Suppose the continuous linear operator C : WN ! Q is defined by ; ðCw qÞ ¼ Cð w; qÞ
q 2 Q: 8 w 2 WN ;
ð34Þ
Then Ker C ¼ SN :
ð35Þ
Þ is coercive in Ker C, i.e., there exists aN > 0 such that Lemma 2. a ð w; w Þ P aN k að w; w wk2WN
8 w 2 Ker C:
ð36Þ
This lemma comes directly from Theorem 1 and Lemma 1. Lemma 3. The bilinear form Cð w; qÞ satisfies the inf–sup condition, namely there exists b > 0 such that sup 2WN nf0g w
jCð w; qÞj P bk qkQ k wkWN
8 q 2 Q n f0g:
ð37Þ
Proof. For a given q 2 Q n f0g, let vk ¼ qkþ1 ; k ¼ 1; 2; . . . ; L 1, vL ¼ 0 and v0 be the unique solution of the following boundary value problem: Dv0 ¼ 0 in Xi ; v0 jCe ¼ q1 ; v0 jCi ¼ 0: Then v ¼ ðv0 ; v1 ; v2 ; . . . ; vL Þ 2 WN n f0g and Cðv; qÞ ¼
Z 0
2p
L X oqk 2 k¼1
oh
þ
L1 X k¼1
! qk qkþ1 dh:
2018
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
For the given quadratic form L1 X
qk qkþ1 6 cL
k¼1
L X
PL1 k¼1
qk qkþ1 , there is a positive number cL < 1 such that
q2k
ð38Þ
k¼1
on the other hand using Parseval equation 2 Z 2p Z 2p oqk q2k dh 6 dh oh 0 0 therefore Cðv; qÞ P
L Z X k¼1
2p
0
oqk oh
cL ðqk Þ
2
ð39Þ
dh P ð1 cL Þ
L Z X k¼1
2p
0
oqk oh
2
2
dh ¼ ð1 cL ÞkqkQ
furthermore 2
2
kvkWN ¼ kv0 k1;Xi þ
L X
2
2
kqk k1;Ce 6 dk qkQ ;
k¼1
where d > 0 is a constant. 8 q 2 Q n f0g, sup 2WN nf0g w
jCð w; qÞj jCðv; qÞj ð1 cL Þ P P k qkQ k wkWN kvkWN d1=2
the inequality follows with b ¼ ð1 cL Þ=d1=2 . By Lemmas 1–3, the following theorem holds. Theorem 2. The mixed variational problem (33) has a unique solution ðuN ; pN Þ 2 WN Q. 4. Finite element approximation of the mixed variational problem (33) Suppose that V1h is a finite element subspace of V1 and S h V1h jCe \ H01 ðCe Þ is a finite element subspace of L L with the mesh size h [6]. Furthermore let WNh ¼ V1h ðS h Þ , Qh ¼ ðS h Þ . Replacing WN and Q in the problem (33) by the two finite element subspaces WNh and Qh gives: Find ð uhN ; pNh Þ 2 WNh Qh such that
H01 ðCe Þ
h Þ þ Cð uhN ; w wh ; pNh Þ ¼ f ð wh Þ a ð Cð uhN ; qh Þ ¼ 0
8 wh 2 WNh ;
8 qh 2 Q h :
ð40Þ
Lemma 4. Let the continuous linear operator Ch : WNh ! Qh be defined by h; qh Þ ¼ Cð wh ; qh Þ ðCh w
8 wh 2 WNh ; qh 2 Q h :
ð41Þ
Then Ker Ch ¼
h
2 w
Z
WNh
0
2p
whk q dh
¼
Z
2p 0
owhk1 oq h dh; k ¼ 1; 2; . . . ; L; 8q 2 S : oh oh
ð42Þ
Lemma 5. The linear operator Ch satisfies the inf–sup condition, namely there exists a constant b > 0 such that
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Cð wh ; qh Þ P bk qh kQ h k w kW N
sup h 2WNh nf0g w
8 qh 2 Q h :
2019
ð43Þ
The proof of Lemma 5 is similar to the proof of Lemma 3, which is omitted here. Lemma 6. For N ¼ 1; 3; 5, i.e., L ¼ 0; 1; 2, the bilinear form a ð; Þ is coercive in Ker Ch , namely there is a positive constant cN such that Þ P cN k a ð w; w wk2WN
8 w 2 Ker Ch :
ð44Þ
Þ, w; w Proof. By the definition of a ð L X
; w Þ ¼ aðw0 ; w0 Þ þ a ðw
aN2k
þ
L X
2p
w2k dh
þ
0
k¼1 2 P akw0 k1;Xi
Z
L X
aN2k
Z
2p
w2k dh
þ
L X
Z
2p
0
k¼1
0
k¼1
aN2kþ1 aN2kþ1
Z
2p
0
k¼1
owk oh owk oh
2 dh 2 dh:
Denote I
L X
aN2k
Z
2p
w2k dh þ 0
k¼1
L X
aN2kþ1
Z
2p
0
k¼1
owk oh
2 dh:
2 Ker Ch , if the following is proved For w 2 L Z 2p X owk dh I P cN oh 0 k¼0
ð45Þ
with cN > 0, the inequality (44) follows directly. 2 Ker Ch , the following holds: For w Z 2p Z 2p owk1 owk dh 8k ¼ 1; 2; . . . ; L; w2k dh ¼ oh oh 0 0 Z
2p
0
owj1 owk1 dh ¼ oh oh
Z
2p 0
ð46Þ
owj2 owk2 dh 81 6 j1 þ k1 ¼ j2 þ k2 : oh oh
ð47Þ
When N ¼ 1, L ¼ 0, obviously the inequality (44) holds. 2 Ker Ch , using (46) When N ¼ 3, L ¼ 1, for w I¼
a32
Z
2p
w21 dh
0
¼
Z 0
2p
a31
þ
ow0 oh
a31 2
Z
2p
0
þ
a33
ow0 oh ow1 oh
2 dh þ 2
a33
Z
2p
0
a32
ow0 ow1 oh oh
ow1 oh !
2 dh
dh;
where fa3k ; k ¼ 1; 2; 3g is given in Table 1. Denote II ¼ a31 x21 þ a33 x22 a32 x1 x2 and then II is a positive definite quadratic form. Take c3 ¼ 3:9E 3 such that II P c3 ðx21 þ x22 Þ, and the inequality (45) follows immediately.
2020
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2 Ker Ch , using (46) and (47) When N ¼ 5, L ¼ 2, for w ! 2 2 2 Z 2p 5 ow0 5 ow1 5 ow2 5 2 5 2 I¼ a1 þ a5 þ a2 w1 þ a4 w2 dh þ a3 oh oh oh 0 ! 2 2 2 Z 2p 5 ow0 5 ow1 5 ow2 5 ow0 ow1 5 ow1 ow2 a4 a1 þ a5 a2 ¼ þ a3 dh oh oh oh oh oh oh oh 0 ! 2 2 2 Z 2p ow ow ow ow ow ow ow ow ow 0 1 2 0 1 1 2 0 2 a54 c a51 þ a55 a52 þ ða53 þ cÞ ¼ dh: oh oh oh oh oh oh oh oh oh 0 Consider the quadratic form II ¼ a51 x21 þ ða53 þ cÞx22 þ a55 x23 a52 x1 x2 a54 x2 x3 cx1 x3
ð48Þ
with c ¼ 0:009579 and then a computation shows that II is a positive definite quadratic form and there is a constant c5 ¼ 2:5E 7 > 0 such that
ð49Þ II P c5 x21 þ x22 þ x23 : The proof is complete.
Theorem 3. For N ¼ 1; 3; 5, the approximate problem (40) has a unique solution ðuhN ; pNh Þ and the following error estimate holds: h h h h kWN þ inf kpN q kQ ; k uN uN kWN þ k pN pN kQ 6 c inf k uN w ð50Þ 2WNh w
qh 2Qh
pN Þ is the solution of mixed variational problem (33) and c is a constant independent of h. where ð uN ; This theorem comes directly from the standard analysis of the mixed finite element method in [5] and Lemmas 1–6.
5. Numerical examples Three examples are given in this section to show the performance of our method. Denote Err as the maximum error on the grid points. 5.1. Example 1 Consider Poisson equation in the exterior domain outside of a circle with the radius 1. The problem is described as the following: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o Du ¼ 0 in X ¼ ðx; yÞ x2 þ y 2 > 1 ; ð51Þ uðx; yÞ ¼
x 0:5 2
ðx 0:5Þ þ y 2
u is bounded at infinity:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o on C ¼ ðx; yÞ x2 þ y 2 ¼ 1 ;
ð52Þ ð53Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2021
The exact solution of this problem is x 0:5 : uðx; yÞ ¼ ðx 0:5Þ2 þ y 2 Introduce a circular artificial boundary Ce with the radius 2. Fig. 1 shows the method of mesh generation. The final mesh is denoted as MeshðNr ; Nh Þ. Fig. 2 shows the relation between the mesh size and the error for different kinds of local artificial boundary conditions N ¼ 1; 3 and 5. 5.2. Example 2 Consider the ideal irrotational incompressible flow around an elliptical obstacle with long axis a ¼ 1 and pffiffiffi short axis b ¼ 0:5. The focus c ¼ 3=2. Suppose the incoming flow is parallel to the long axis. Denote the
Fig. 1. Mesh plot outside of a circle. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 24.
Fig. 2. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Let Nh ¼ 6Nr in the implementation.
2022
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
surface of the obstacle as C and the exterior domain outside of C as X, and the problem can be described with the velocity field ~ v ¼ ðu; vÞ as the following: r ~ v ¼ 0 in X; r ~ v ¼ 0 in X; ~ v n ¼ 0 on C; Z ~ v d~ l ¼ c; C
~ v ! ðV ; 0Þ at infinity: Here, c denotes the circulation of flow around the obstacle and n is the exterior normal direction. Because the flow is incompressible, the stream function w can be induced, which is related with ~ v by ow ¼ u; oy
ow ¼ v: ox
Then w satisfies the following equations: Dw ¼ 0 in X; w¼0 w
on C;
c ln r Vy is bounded at infinity: 2p
The ideal irrotational flow can be considered as a combination of a flow without circulation and a flow with zero velocity at infinity. So two cases, c ¼ 2p, V ¼ 0 and c ¼ 0, V ¼ 1, will be considered in the following. It is known that for c ¼ 2p, V ¼ 0, incompressible flow around an ellipse with zero the stream function is given by
Fig. 3. Mesh plot outside of an ellipse. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 24.
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2023
Fig. 4. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Suppose the circulation is 2p and the velocity at infinity is zero.
a W w ¼ Real cosh1 cosh1 c c and for c ¼ 0, V ¼ 1, the stream function is given by 1 1 1 c Imagine ecosh ðW =cÞ þ e2 cosh ða=cÞcosh ðW =cÞ ; w¼ 2 where W ¼ x þ iy: Introduce a circular artificial boundary with radius 2. Fig. 3 shows the mesh generation. The final mesh is denoted as MeshðNr ; Nh Þ. Figs. 4 and 6 show the relation between the mesh size and the error for different kinds of local artificial boundary conditions N ¼ 1; 3 and 5. Figs. 5 and 7 show the stream line.
Fig. 5. The stream line plot outside of an ellipse with velocity zero at infinity and circulation 2p. In the implementation, Nr ¼ 64, Nh ¼ 384, N ¼ 5.
2024
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Fig. 6. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Suppose the circulation is zero and the velocity at infinity is one.
Fig. 7. The stream line plot outside of an ellipse with circulation zero and velocity one at infinity. In the implementation, Nr ¼ 64, Nh ¼ 384, N ¼ 5.
5.3. Example 3 Now consider the ideal irrotational incompressible flow around a 2 2 rectangular obstacle. Let the origin be the center of the square. Introduce a circular artificial boundary with radius 4. Fig. 8 shows the method of mesh generation. Fig. 9 shows the stream line when the velocity of flow is zero at infinity and Fig. 10 shows the stream line when the circulation is zero. In both cases, Nr ¼ 32, Nh ¼ 248 and local artificial boundary condition N ¼ 5 are used in the computation. 5.4. Experimental observations The following can be observed in Figs. 2, 4 and 6: • the mixed finite element method is effective to deal with the higher-order local artificial boundary conditions;
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2025
Fig. 8. Mesh plot outside of a square. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 32.
• higher-order local artificial boundary conditions show their advantages over the lower-order ones when the mesh is properly refined; • for each local artificial boundary condition, the error remains nearly unchangeable after the mesh is sufficiently refined. The error in the numerical implementation has originated from two sources: one is the approximation of the artificial boundary conditions, and the other is the approximation of finite element method. When the low-order local artificial boundary condition is used, the main part of the error comes from the approximation of the artificial boundary condition. Then a more refined mesh cannot reduce the error efficiently any more; in order to obtain a numerical solution with a much higher accuracy, higher-order local artificial boundary conditions are necessary. The numerical results are compatible with the analysis.
Fig. 9. The steam line plot outside of a square with circulation 2p and velocity zero at infinity. In the computation, Nr ¼ 32, Nh ¼ 248, N ¼ 5.
2026
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Fig. 10. The stream line plot outside of a square with circulation zero and velocity one at infinity. In the computation, Nr ¼ 32, Nh ¼ 248, N ¼ 5.
6. Conclusions When high-order local artificial boundary conditions are used, the new unknowns on artificial boundary are introduced. Then the problem with high-order derivatives of unknown function on the artificial boundary is reduced to an equivalent mixed variational problem in which there are no high-order derivatives involved. It has been proved that for some high-order local artificial boundary conditions, the discrete mixed variational formulations are well-posed. The error estimate is also obtained. This method avoids the construction of high-order finite elements on the adjacent layer of the artificial boundary, which seems to be very complicated when high-order local artificial boundary conditions are used. Two remarks can be made up here on local artificial boundary conditions: • when N > 0 is an even number, the bilinear form bN ð; Þ loses coercivity and the stability of the reduced problem cannot be guaranteed, and so these local artificial boundary conditions are suggested as not to be used; • the local artificial boundary condition when N ¼ 5 is accurate enough for the practical use. High-order local artificial boundary condition can bring up an ill-conditioned stiffness matrix.
References [1] R.A. Adams, Sobolev Spaces, 1975. [2] W. Bao, Artificial boundary conditions for two-dimensional incompressible viscous flows around an obstacle, Comput. Methods Appl. Mech. Engrg. 147 (1997) 263. [3] W. Bao, H. Han, Local artificial boundary conditions for the incompressible viscous flow in a slip channel, J. Comput. Math. 15 (1997) 335. [4] W. Bao, H. Han, High-order local artificial boundary conditions for problems in unbounded domains, Comput. Methods Appl. Mech. Engrg. 188 (2000) 455–471. [5] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [6] P.G. Ciarlet, the Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [7] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (1977) 629–651. [8] K. Feng, Asymptotic radiation conditions for reduced wave equations, J. Comput. Math. 2 (1984) 130–138.
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2027
[9] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992. [10] D. Givoli, J.B. Keller, Special finite elements for use with high-order boundary conditions, Comput. Methods Appl. Mech. Engrg. 119 (1994) 199–213. [11] D. Givoli, I. Patlashenko, An optimal high-order non-reflecting finite element scheme for wave scattering problems, Comput. Methods Appl. Mech. Engrg., to appear. [12] D. Givoli, I. Patlashenko, J. Keller, High-order boundary conditions and finite elements for infinite domains, Comput. Methods Appl. Mech. Engrg. 143 (1997) 13–39. [13] D. Givoli, I. Patlashenko, J. Keller, Discrete Dirichlet-to-Neumann maps for unbounded domains, Comput. Methods Appl. Mech. Engrg. 164 (1998) 173–185. [14] C.I. Goldstein, A finite element method for solving Helmholtz type equations in waveguides and other unbounded domains, Math. Comp. 39 (1982) 309–324. [15] T.M. Hagstrom, H.B. Keller, Exact boundary conditions at artificial boundary for partial differential equations in cylinders, SIAM J. Math. Anal. 17 (1986) 322–341. [16] T.M. Hagstrom, H.B. Keller, Asymptotic boundary conditions and numerical methods for nonlinear elliptic problems on unbounded domains, Math. Comp. 48 (1987) 449–470. [17] L. Halpern, M. Schatzman, Artificial boundary conditions for incompressible viscous flows, SIAM J. Math. Anal. 20 (1989) 308– 353. [18] H. Han, W. Bao, An artificial boundary condition for the incompressible viscous flows in a no-slip channel, J. Comput. Math. 13 (1995) 51–63. [19] H. Han, W. Bao, The artificial boundary conditions for incompressible viscous flows in channels, Lecture Notes Numer. Appl. Anal. 14 (1995) 47. [20] H. Han, W. Bao, An artificial boundary condition for two-dimensional incompressible viscous flows using the method of lines, Int. J. Numer. Methods Fluids 22 (1996) 483–493. [21] H. Han, J. Lu, W. Bao, A discrete artificial boundary condition for steady incompressible viscous flows in a no-slip channel using a fast iterative method, J. Comput. Phys. 114 (1994) 201–208. [22] H. Han, W. Bao, T. Wang, Numerical simulation for the problem of infinite elastic foundation, Comput. Methods Appl. Mech. Engrg. 147 (1997) 369–385. [23] H. Han, X. Wu, Approximation of infinite boundary condition and its application to finite element method, J. Comput. Math. 3 (1985) 179–192. [24] H. Han, X. Wu, The approximation of exact boundary condition at an artificial boundary for linear elastic equation and its application, Math. Comp. 59 (1992) 21–27. [25] J.L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Springer, Berlin, 1971. [26] D. Yu, Approximation of boundary conditions at infinity for a harmonic equation, J. Comput. Math. 3 (1985) 220–227.
Mixed finite element method and high-order local artificial boundary conditions for exterior problems of elliptic equation q Houde Han *, Chunxiong Zheng Department of Mathematical Sciences, Tsinghua University, Beijing 100084, People’s Republic of China Received 21 February 2001; received in revised form 18 October 2001; accepted 19 October 2001
Abstract The mixed finite element method is used to solve the exterior elliptic problem with high-order local artificial boundary conditions. New unknowns are introduced to reduce the order of the derivatives to two. This leads to an equivalent mixed variational problem such that the normal finite element can be used and special finite elements are no longer needed on the adjacent layer of the artificial boundary. Error estimates are obtained for some local artificial boundary conditions with prescribed order. Numerical examples are presented and the results demonstrate the effectiveness of this method. 2001 Elsevier Science B.V. All rights reserved. Keywords: Mixed finite element method; Artificial boundary condition; Elliptic equation
1. Introduction Elliptic problems on exterior domain arise in many application fields. For example, the incompressible irrotational flow around a body is described by the exterior problems of Laplace equation. There are many numerical methods to solve these exterior problems. The most popular is the artificial boundary method. In the procedure of this method, an artificial boundary is firstly introduced to divide the exterior domain into a bounded part and an unbounded part. After a suitable artificial boundary condition is set up on the artificial boundary, the original exterior problem is reduced to a new one only defined on the bounded domain, and then it can be solved by a proper numerical method. In the last two decades, many authors have worked on this subject for various problems by different techniques (see [7–9,12–25] and the references
q
This work was supported partly by the Special Funds for Major State Basic Research Projects of China and the National Natural Science Foundation of China. Computation was supported by the State Key Laboratory of the Scientific and Engineering Computing in China. * Corresponding author. E-mail address: [email protected] (H. Han). 0045-7825/01/$ - see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 3 6 5 - 6
2012
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
therein). Generally speaking, the artificial boundary conditions can be classified into nonlocal artificial boundary conditions, local artificial boundary conditions and discrete artificial boundary conditions. A thorough review concerning nonlocal and local artificial boundary conditions can be found in [12]. For the discrete artificial boundary conditions, the readers can refer to [2,3,13,18–22]. For high-order local artificial boundary conditions, since the appearance of high-order derivatives, the new difficulty is posed in the implementation with a conforming finite element method. To allow the use of high-order local artificial boundary conditions, Givoli and others [10,12] developed a family of special elements and in [11], Givoli and Patlashenko proposed a new finite element scheme for the numerical simulation of time-harmonic wave scattering problems in unbounded domains. In this paper, mixed finite element method will be used to overcome this new difficulty.
2. High-order local artificial boundary conditions Suppose Ci is a bounded, simple closed curve in R2 and X is the exterior domain with boundary Ci . The boundary Ci is decomposed into two disjoint parts: CD with a nonzero measure where a Dirichlet boundary condition is given and CN where a Neumann boundary condition is given. Consider the following exterior problems of a linear elliptic second-order equation: r ðjðxÞruðxÞÞ þ bðxÞuðxÞ ¼ f ðxÞ in X;
ð1Þ
u¼g
ð2Þ
on CD ;
ou ¼k on
ð3Þ
on CN ;
u is bounded
when jxj ! 1;
ð4Þ
, g and k are the where uðxÞ is the unknown function, jðxÞ P 1, bðxÞ P 0, f ðxÞ are the given functions on X given functions on CD and CN . Furthermore, suppose suppðf Þ, suppðj 1Þ and suppðbÞ are compact, jðxÞ 2 L1 ðXÞ, bðxÞ 2 L1 ðXÞ, f ðxÞ 2 L2 ðXÞ, g 2 H 3=2 ðCD Þ, kðxÞ 2 H 1=2 ðCN Þ. For brevity let gðxÞ 0 in the following. Let R0 ¼ 0 þ maxfr: ðr; hÞ 2 Ci [ suppðf Þ [ suppðjÞ [ suppðbÞg; where 0 is a positive constant and C0 ¼ fðR0 ; hÞ: 0 6 h 6 2pg. Then C0 X. Introduce a circular artificial boundary Ce fðR; hÞ: 0 6 h 6 2pg with R > R0 . Ce divides the exterior domain into two parts: the e . On the unbounded unbounded part Xe ¼ fðr; hÞ: r > R; 0 6 h 6 2pg and the bounded part Xi ¼ X n X domain Xe , jðxÞ 1, bðxÞ 0, f ðxÞ 0 and Eq. (1) is reduced to a Laplace equation. Now recall the highorder local artificial boundary conditions on the given artificial boundary Ce [4,12,23]. Suppose uðxÞ is the solution of problem (1)–(4). Then on X0 fðr; hÞ: r > R0 ; 0 6 h 6 2pg, uðxÞ satisfies Laplace equation and condition (4). By separation of variables, uðxÞ can be given on X0 by n 1 a0 X R0 ðan cos nh þ bn sin nhÞ ð5Þ uðr; hÞ ¼ þ 2 r n¼1 with 1 an ¼ p
Z 0
2p
1 uðR0 ; hÞ cos nh dh ¼ p
R R0
n Z
2p
uðR; hÞ cos nh dh; 0
ð6Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
bn ¼
1 p
Z
2p
uðR0 ; hÞ sin nh dh ¼ 0
1 p
R R0
n Z
2013
2p
ð7Þ
uðR; hÞ sin nh dh: 0
Therefore n 1 ou ouðR; hÞ 1X R0 ¼ n ðan cos nh þ bn sin nhÞ on Ce or R n¼1 R Z 2p 1 1 X n uðR; h0 Þ cos nðh h0 Þ dh0 : ¼ pR n¼1 0
ð8Þ
This is the desired exact boundary condition on Ce . Furthermore, n 1 a0 X R0 ðan cos nh þ bn sin nhÞ uðR; hÞ ¼ þ 2 R n¼1
ð9Þ
and n 1 o2m uðR; hÞ X 2 m R0 ¼ ðn Þ ðan cos nh þ bn sin nhÞ R oh2m n¼1
ð10Þ
with positive integer m. A computation shows
N N 2m ouðR; hÞ 1 X o2m uðR; hÞ ouðR; hÞ 1 X m m N o uðR; hÞ ð1Þ aNm ¼ ð1Þ a m on R m¼1 or R m¼1 oh2m oh2m n 1 1X R0 ¼ An ðan cos nh þ bn sin nhÞ; R n¼1 R
where N is a positive integer and faN1 ; aN2 ; . . . ; aNN g are N constants and An ¼ n
N X
aNm n2m ;
n ¼ 1; 2; . . .
m¼1
Let faN1 ; aN2 ; . . . ; aNN g be the solution of the following linear algebraic equations: N X
aNm n2m ¼ n;
n ¼ 1; 2; . . . ; N :
ð11Þ
m¼1
Then
n N 1 ouðR; hÞ 1 X o2m uðR; hÞ 1 X R0 1 m ð1Þ aNm ¼ A ða cos nh þ b sin nhÞ ¼ O : n n n on R m¼1 R n¼1 RN þ2 R oh2m
Table 1 Solution of algebraic equations (11) aN1 N ¼1 N ¼2 N ¼3 N ¼4 N ¼5
aN2
aN3
aN4
1 60 11 360 71 1728
1 1008
aN5
1 7 6 37 30 533 420 1627 1260
16 14 43 144 107 324
13 6048
1 25 920
ð12Þ
2014
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Omitting the RHS of (12) gives the following high-order local artificial boundary conditions on the given artificial boundary, which are exact for the first N modes: N 2m ou 1X m N o uðR; hÞ ¼ ð1Þ a : ð13Þ m on Ce R m¼1 oh2m For any positive integer N, it is easy to be verified that algebraic equations (11) have a unique solution. The solution for N ¼ 1; 2; . . . ; 5 is illustrated in Table 1. Applying the local artificial boundary conditions gives a series of approximate boundary value problems defined on the bounded domain Xi : r ðjðxÞruN ðxÞÞ þ bðxÞuN ðxÞ ¼ f ðxÞ
ð14Þ
in Xi ;
uN ¼ 0 on CD ;
ð15Þ
ouN ¼k on
ð16Þ
on CN ;
N ouN 1X o2m uN m ¼ ð1Þ aNm R m¼1 on oh2m
ð17Þ
on Ce :
Let H m ðXi Þ and H s ðCe Þ be the usual Sobolev spaces on Xi and Ce with integer m and real number s [1]. Suppose H1 ðXi Þ ¼ fv 2 H 1 ðXi Þ: v ¼ 0 on CD g: Then the weak form of problem (1)–(3) with the exact boundary condition (8) is: Find u 2 H1 ðXi Þ such that aðu; vÞ þ bðu; vÞ ¼ f ðvÞ 8v 2 H1 ðXi Þ; where aðu; vÞ ¼
Z
jðxÞruðxÞ rvðxÞ dx þ
Xi
Z
ð18Þ
bðxÞuðxÞvðxÞ dx 8u; v 2 H 1 ðXi Þ;
Z Z 1 X n 2p 2p cos nðh h0 ÞuðR; hÞvðR; h0 Þ dh dh0 bðu; vÞ ¼ p 0 0 n¼1 f ðvÞ ¼
Z
ð19Þ
Xi
f ðxÞvðxÞ dx þ
Xi
Z
8u; v 2 H 1 ðXi Þ;
jðxÞkðxÞvðxÞ ds 8v 2 H 1 ðXi Þ:
ð20Þ
ð21Þ
CN
In order to formulate the weak form of the problem (14)–(17), define VN ¼ fv 2 H1 ðXi Þ: vjCe 2 H N ðCe Þg with norm kvkVN jvj1;Xi þ jvjN ;Ce :
ð22Þ
Let bN ðu; vÞ ¼
N X m¼1
aNm
Z
2p 0
om uðR; hÞ om vðR; hÞ dh: ohm ohm
ð23Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2015
Then the boundary value problem (14)–(17) is equivalent to the following variational problem: Find uN 2 VN such that aðuN ; vÞ þ bN ðuN ; vÞ ¼ f ðvÞ
8v 2 VN :
ð24Þ
Theorem 1. Suppose N is an odd number: • The bilinear form að; Þ þ bN ð; Þ is coercive on VN and the variational problem (24) is well-posed. • Let uðxÞ be the solution of problem (1)–(4). Then the restriction of uðxÞ on Xi is the solution of variational problem (18). Assume that uN is the solution of the problem (24). Then the following error estimate holds: ku uN kVN 6 CN
R0 R
N þ1 jujN ;C0 ;
ð25Þ
where CN is a constant only depending on N.
Remark. In [4], this theorem is proved under the hypothesis of 1 6 N < 20. But the restriction on N is not essential. The key point is to prove the coercivity property of bN ð; Þ on H N ðCe Þ, which has been done in [26] for all odd numbers. So the proof of this theorem is omitted here. When N > 1, the high-order derivatives of the unknown function on boundary Ce appear in the variational problem (24), therefore, the standard C ð0Þ -element cannot be used to get the numerical approximation of problem (24). Givoli and others [10,12] developed a family of special elements and their extension, a hierarchy of finite elements, to overcome this new difficulty. In the following sections, a different approach, mixed finite element method, will be proposed to compute the numerical solution of the boundary value problem (14)–(17) for N > 1.
3. Mixed variational formulation of problem (14)–(17) Suppose that N ¼ 2L þ 1 is an odd number and J ðwÞ ¼ 12ðaðw; wÞ þ bN ðw; wÞÞ f ðwÞ
8w 2 VN :
According to Theorem 1, the boundary value problem (14)–(17) is equivalent to the following minimization problem: Find uN 2 VN such that J ðuN Þ ¼ inf J ðwÞ:
ð26Þ
w2VN
For 8w 2 VN , introduce new unknowns on the artificial boundary Ce denoted by wk ðR; hÞ ¼
o2k wðR; hÞ ; oh2k
k ¼ 1; 2; . . . ; L
ð27Þ
and w0 ðr; hÞ ¼ wðr; hÞ
8ðr; hÞ 2 Xi :
ð28Þ
2016
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
The bilinear form bN ðw; wÞ can be rewritten as follows: 2 Z 2p m N X o wðR; hÞ N bN ðw; wÞ ¼ am dh ohm 0 m¼1 2 2 Z 2p 2kþ1 Z 2p 2k L L X X o wðR; hÞ o wðR; hÞ N N ¼ a2kþ1 dh þ a2k dh oh2kþ1 oh2k 0 0 k¼0 k¼1 2 Z 2p Z 2p L L X X owk ðR; hÞ N ðwk ðR; hÞÞ2 dh: ¼ a2kþ1 dh þ aN2k oh 0 0 k¼0 k¼1 Denote H01 ðCe Þ ¼
w 2 H 1 ðCe Þ and
Z
w¼0 ;
Ce
n o ¼ ðw0 ; w1 ; w2 ; . . . ; wL Þ j w0 2 V1 ; wk 2 H01 ðCe Þ; k ¼ 1; 2; . . . ; L ; WN ¼ w n o ¼ ðq1 ; q2 ; . . . ; qL Þ j qk 2 H01 ðCe Þ; k ¼ 1; 2; . . . ; L Q¼ q and Þ ¼ bN ð u; w
L X
aN2kþ1
Z 0
k¼0
2p
L X ouk ðR; hÞ owk ðR; hÞ dh þ aN2k oh oh k¼1
Z
2p
uk ðR; hÞwk ðR; hÞ dh;
0
Þ ¼ aðu0 ; w0 Þ þ bN ð Þ; a ð u; w u; w f ð wÞ ¼ f ðw0 Þ; Þ f ð J ð wÞ ¼ 12a ð w; w wÞ; Cð w; qÞ ¼
L Z X k¼1
2p 0
owk1 oqk w k qk þ dh oh oh
8 w 2 WN ; 8q 2 Q;
2 WN and q 2 Q. Then where u 2 WN , w Þ w; w bN ðw; wÞ ¼ bN ð and Eq. (27) are equivalent to the following variational form Cð w; qÞ ¼ 0
8 q 2 Q:
ð29Þ
Furthermore introduce n o jw 2 WN and Cð w; qÞ ¼ 0; 8 q2Q : SN ¼ w Obviously SN is a subspace of WN and the minimization problem (26) is equivalent to the following conditional minimization problem: Find u 2 SN such that uÞ ¼ min J ð wÞ: J ð 2SN w
ð30Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2017
Introduce the Lagrange multiplier q 2 Q and let Þ þ Cð Ið w; qÞ ¼ 12a ð w; w w; qÞ f ð wÞ:
ð31Þ
Then the conditional minimization problem (30) is equivalent to the following saddle problem [5]: Find ð uN ; pN Þ 2 WN Q such that pN Þ ¼ inf sup Ið w; qÞ: Ið uN ;
ð32Þ
2WN w q2Q
Furthermore the saddle problem is equivalent to the following mixed variational problem: Find ð uN ; pN Þ 2 WN Q such that Þ þ Cð a ð uN ; w w; pN Þ ¼ f ð wÞ 8 w 2 WN ; Þ ¼ 0 8 Cð uN ; q q 2 Q:
ð33Þ
Finally the mixed variational problem (33) is obtained, which is equivalent to the boundary value problem (14)–(17). And now there are no high-order derivatives of the unknowns involved in the formulation. Lemma 1. Suppose the continuous linear operator C : WN ! Q is defined by ; ðCw qÞ ¼ Cð w; qÞ
q 2 Q: 8 w 2 WN ;
ð34Þ
Then Ker C ¼ SN :
ð35Þ
Þ is coercive in Ker C, i.e., there exists aN > 0 such that Lemma 2. a ð w; w Þ P aN k að w; w wk2WN
8 w 2 Ker C:
ð36Þ
This lemma comes directly from Theorem 1 and Lemma 1. Lemma 3. The bilinear form Cð w; qÞ satisfies the inf–sup condition, namely there exists b > 0 such that sup 2WN nf0g w
jCð w; qÞj P bk qkQ k wkWN
8 q 2 Q n f0g:
ð37Þ
Proof. For a given q 2 Q n f0g, let vk ¼ qkþ1 ; k ¼ 1; 2; . . . ; L 1, vL ¼ 0 and v0 be the unique solution of the following boundary value problem: Dv0 ¼ 0 in Xi ; v0 jCe ¼ q1 ; v0 jCi ¼ 0: Then v ¼ ðv0 ; v1 ; v2 ; . . . ; vL Þ 2 WN n f0g and Cðv; qÞ ¼
Z 0
2p
L X oqk 2 k¼1
oh
þ
L1 X k¼1
! qk qkþ1 dh:
2018
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
For the given quadratic form L1 X
qk qkþ1 6 cL
k¼1
L X
PL1 k¼1
qk qkþ1 , there is a positive number cL < 1 such that
q2k
ð38Þ
k¼1
on the other hand using Parseval equation 2 Z 2p Z 2p oqk q2k dh 6 dh oh 0 0 therefore Cðv; qÞ P
L Z X k¼1
2p
0
oqk oh
cL ðqk Þ
2
ð39Þ
dh P ð1 cL Þ
L Z X k¼1
2p
0
oqk oh
2
2
dh ¼ ð1 cL ÞkqkQ
furthermore 2
2
kvkWN ¼ kv0 k1;Xi þ
L X
2
2
kqk k1;Ce 6 dk qkQ ;
k¼1
where d > 0 is a constant. 8 q 2 Q n f0g, sup 2WN nf0g w
jCð w; qÞj jCðv; qÞj ð1 cL Þ P P k qkQ k wkWN kvkWN d1=2
the inequality follows with b ¼ ð1 cL Þ=d1=2 . By Lemmas 1–3, the following theorem holds. Theorem 2. The mixed variational problem (33) has a unique solution ðuN ; pN Þ 2 WN Q. 4. Finite element approximation of the mixed variational problem (33) Suppose that V1h is a finite element subspace of V1 and S h V1h jCe \ H01 ðCe Þ is a finite element subspace of L L with the mesh size h [6]. Furthermore let WNh ¼ V1h ðS h Þ , Qh ¼ ðS h Þ . Replacing WN and Q in the problem (33) by the two finite element subspaces WNh and Qh gives: Find ð uhN ; pNh Þ 2 WNh Qh such that
H01 ðCe Þ
h Þ þ Cð uhN ; w wh ; pNh Þ ¼ f ð wh Þ a ð Cð uhN ; qh Þ ¼ 0
8 wh 2 WNh ;
8 qh 2 Q h :
ð40Þ
Lemma 4. Let the continuous linear operator Ch : WNh ! Qh be defined by h; qh Þ ¼ Cð wh ; qh Þ ðCh w
8 wh 2 WNh ; qh 2 Q h :
ð41Þ
Then Ker Ch ¼
h
2 w
Z
WNh
0
2p
whk q dh
¼
Z
2p 0
owhk1 oq h dh; k ¼ 1; 2; . . . ; L; 8q 2 S : oh oh
ð42Þ
Lemma 5. The linear operator Ch satisfies the inf–sup condition, namely there exists a constant b > 0 such that
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Cð wh ; qh Þ P bk qh kQ h k w kW N
sup h 2WNh nf0g w
8 qh 2 Q h :
2019
ð43Þ
The proof of Lemma 5 is similar to the proof of Lemma 3, which is omitted here. Lemma 6. For N ¼ 1; 3; 5, i.e., L ¼ 0; 1; 2, the bilinear form a ð; Þ is coercive in Ker Ch , namely there is a positive constant cN such that Þ P cN k a ð w; w wk2WN
8 w 2 Ker Ch :
ð44Þ
Þ, w; w Proof. By the definition of a ð L X
; w Þ ¼ aðw0 ; w0 Þ þ a ðw
aN2k
þ
L X
2p
w2k dh
þ
0
k¼1 2 P akw0 k1;Xi
Z
L X
aN2k
Z
2p
w2k dh
þ
L X
Z
2p
0
k¼1
0
k¼1
aN2kþ1 aN2kþ1
Z
2p
0
k¼1
owk oh owk oh
2 dh 2 dh:
Denote I
L X
aN2k
Z
2p
w2k dh þ 0
k¼1
L X
aN2kþ1
Z
2p
0
k¼1
owk oh
2 dh:
2 Ker Ch , if the following is proved For w 2 L Z 2p X owk dh I P cN oh 0 k¼0
ð45Þ
with cN > 0, the inequality (44) follows directly. 2 Ker Ch , the following holds: For w Z 2p Z 2p owk1 owk dh 8k ¼ 1; 2; . . . ; L; w2k dh ¼ oh oh 0 0 Z
2p
0
owj1 owk1 dh ¼ oh oh
Z
2p 0
ð46Þ
owj2 owk2 dh 81 6 j1 þ k1 ¼ j2 þ k2 : oh oh
ð47Þ
When N ¼ 1, L ¼ 0, obviously the inequality (44) holds. 2 Ker Ch , using (46) When N ¼ 3, L ¼ 1, for w I¼
a32
Z
2p
w21 dh
0
¼
Z 0
2p
a31
þ
ow0 oh
a31 2
Z
2p
0
þ
a33
ow0 oh ow1 oh
2 dh þ 2
a33
Z
2p
0
a32
ow0 ow1 oh oh
ow1 oh !
2 dh
dh;
where fa3k ; k ¼ 1; 2; 3g is given in Table 1. Denote II ¼ a31 x21 þ a33 x22 a32 x1 x2 and then II is a positive definite quadratic form. Take c3 ¼ 3:9E 3 such that II P c3 ðx21 þ x22 Þ, and the inequality (45) follows immediately.
2020
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2 Ker Ch , using (46) and (47) When N ¼ 5, L ¼ 2, for w ! 2 2 2 Z 2p 5 ow0 5 ow1 5 ow2 5 2 5 2 I¼ a1 þ a5 þ a2 w1 þ a4 w2 dh þ a3 oh oh oh 0 ! 2 2 2 Z 2p 5 ow0 5 ow1 5 ow2 5 ow0 ow1 5 ow1 ow2 a4 a1 þ a5 a2 ¼ þ a3 dh oh oh oh oh oh oh oh 0 ! 2 2 2 Z 2p ow ow ow ow ow ow ow ow ow 0 1 2 0 1 1 2 0 2 a54 c a51 þ a55 a52 þ ða53 þ cÞ ¼ dh: oh oh oh oh oh oh oh oh oh 0 Consider the quadratic form II ¼ a51 x21 þ ða53 þ cÞx22 þ a55 x23 a52 x1 x2 a54 x2 x3 cx1 x3
ð48Þ
with c ¼ 0:009579 and then a computation shows that II is a positive definite quadratic form and there is a constant c5 ¼ 2:5E 7 > 0 such that
ð49Þ II P c5 x21 þ x22 þ x23 : The proof is complete.
Theorem 3. For N ¼ 1; 3; 5, the approximate problem (40) has a unique solution ðuhN ; pNh Þ and the following error estimate holds: h h h h kWN þ inf kpN q kQ ; k uN uN kWN þ k pN pN kQ 6 c inf k uN w ð50Þ 2WNh w
qh 2Qh
pN Þ is the solution of mixed variational problem (33) and c is a constant independent of h. where ð uN ; This theorem comes directly from the standard analysis of the mixed finite element method in [5] and Lemmas 1–6.
5. Numerical examples Three examples are given in this section to show the performance of our method. Denote Err as the maximum error on the grid points. 5.1. Example 1 Consider Poisson equation in the exterior domain outside of a circle with the radius 1. The problem is described as the following: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o Du ¼ 0 in X ¼ ðx; yÞ x2 þ y 2 > 1 ; ð51Þ uðx; yÞ ¼
x 0:5 2
ðx 0:5Þ þ y 2
u is bounded at infinity:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n o on C ¼ ðx; yÞ x2 þ y 2 ¼ 1 ;
ð52Þ ð53Þ
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2021
The exact solution of this problem is x 0:5 : uðx; yÞ ¼ ðx 0:5Þ2 þ y 2 Introduce a circular artificial boundary Ce with the radius 2. Fig. 1 shows the method of mesh generation. The final mesh is denoted as MeshðNr ; Nh Þ. Fig. 2 shows the relation between the mesh size and the error for different kinds of local artificial boundary conditions N ¼ 1; 3 and 5. 5.2. Example 2 Consider the ideal irrotational incompressible flow around an elliptical obstacle with long axis a ¼ 1 and pffiffiffi short axis b ¼ 0:5. The focus c ¼ 3=2. Suppose the incoming flow is parallel to the long axis. Denote the
Fig. 1. Mesh plot outside of a circle. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 24.
Fig. 2. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Let Nh ¼ 6Nr in the implementation.
2022
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
surface of the obstacle as C and the exterior domain outside of C as X, and the problem can be described with the velocity field ~ v ¼ ðu; vÞ as the following: r ~ v ¼ 0 in X; r ~ v ¼ 0 in X; ~ v n ¼ 0 on C; Z ~ v d~ l ¼ c; C
~ v ! ðV ; 0Þ at infinity: Here, c denotes the circulation of flow around the obstacle and n is the exterior normal direction. Because the flow is incompressible, the stream function w can be induced, which is related with ~ v by ow ¼ u; oy
ow ¼ v: ox
Then w satisfies the following equations: Dw ¼ 0 in X; w¼0 w
on C;
c ln r Vy is bounded at infinity: 2p
The ideal irrotational flow can be considered as a combination of a flow without circulation and a flow with zero velocity at infinity. So two cases, c ¼ 2p, V ¼ 0 and c ¼ 0, V ¼ 1, will be considered in the following. It is known that for c ¼ 2p, V ¼ 0, incompressible flow around an ellipse with zero the stream function is given by
Fig. 3. Mesh plot outside of an ellipse. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 24.
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2023
Fig. 4. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Suppose the circulation is 2p and the velocity at infinity is zero.
a W w ¼ Real cosh1 cosh1 c c and for c ¼ 0, V ¼ 1, the stream function is given by 1 1 1 c Imagine ecosh ðW =cÞ þ e2 cosh ða=cÞcosh ðW =cÞ ; w¼ 2 where W ¼ x þ iy: Introduce a circular artificial boundary with radius 2. Fig. 3 shows the mesh generation. The final mesh is denoted as MeshðNr ; Nh Þ. Figs. 4 and 6 show the relation between the mesh size and the error for different kinds of local artificial boundary conditions N ¼ 1; 3 and 5. Figs. 5 and 7 show the stream line.
Fig. 5. The stream line plot outside of an ellipse with velocity zero at infinity and circulation 2p. In the implementation, Nr ¼ 64, Nh ¼ 384, N ¼ 5.
2024
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Fig. 6. The change plot of maximum grid error on the mesh size for local artificial boundary condition N ¼ 1, 3 and 5. Suppose the circulation is zero and the velocity at infinity is one.
Fig. 7. The stream line plot outside of an ellipse with circulation zero and velocity one at infinity. In the implementation, Nr ¼ 64, Nh ¼ 384, N ¼ 5.
5.3. Example 3 Now consider the ideal irrotational incompressible flow around a 2 2 rectangular obstacle. Let the origin be the center of the square. Introduce a circular artificial boundary with radius 4. Fig. 8 shows the method of mesh generation. Fig. 9 shows the stream line when the velocity of flow is zero at infinity and Fig. 10 shows the stream line when the circulation is zero. In both cases, Nr ¼ 32, Nh ¼ 248 and local artificial boundary condition N ¼ 5 are used in the computation. 5.4. Experimental observations The following can be observed in Figs. 2, 4 and 6: • the mixed finite element method is effective to deal with the higher-order local artificial boundary conditions;
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2025
Fig. 8. Mesh plot outside of a square. Nr and Nh denote the numbers of divisions in the radial and azimuthal directions. In this plot, Nr ¼ 4, Nh ¼ 32.
• higher-order local artificial boundary conditions show their advantages over the lower-order ones when the mesh is properly refined; • for each local artificial boundary condition, the error remains nearly unchangeable after the mesh is sufficiently refined. The error in the numerical implementation has originated from two sources: one is the approximation of the artificial boundary conditions, and the other is the approximation of finite element method. When the low-order local artificial boundary condition is used, the main part of the error comes from the approximation of the artificial boundary condition. Then a more refined mesh cannot reduce the error efficiently any more; in order to obtain a numerical solution with a much higher accuracy, higher-order local artificial boundary conditions are necessary. The numerical results are compatible with the analysis.
Fig. 9. The steam line plot outside of a square with circulation 2p and velocity zero at infinity. In the computation, Nr ¼ 32, Nh ¼ 248, N ¼ 5.
2026
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
Fig. 10. The stream line plot outside of a square with circulation zero and velocity one at infinity. In the computation, Nr ¼ 32, Nh ¼ 248, N ¼ 5.
6. Conclusions When high-order local artificial boundary conditions are used, the new unknowns on artificial boundary are introduced. Then the problem with high-order derivatives of unknown function on the artificial boundary is reduced to an equivalent mixed variational problem in which there are no high-order derivatives involved. It has been proved that for some high-order local artificial boundary conditions, the discrete mixed variational formulations are well-posed. The error estimate is also obtained. This method avoids the construction of high-order finite elements on the adjacent layer of the artificial boundary, which seems to be very complicated when high-order local artificial boundary conditions are used. Two remarks can be made up here on local artificial boundary conditions: • when N > 0 is an even number, the bilinear form bN ð; Þ loses coercivity and the stability of the reduced problem cannot be guaranteed, and so these local artificial boundary conditions are suggested as not to be used; • the local artificial boundary condition when N ¼ 5 is accurate enough for the practical use. High-order local artificial boundary condition can bring up an ill-conditioned stiffness matrix.
References [1] R.A. Adams, Sobolev Spaces, 1975. [2] W. Bao, Artificial boundary conditions for two-dimensional incompressible viscous flows around an obstacle, Comput. Methods Appl. Mech. Engrg. 147 (1997) 263. [3] W. Bao, H. Han, Local artificial boundary conditions for the incompressible viscous flow in a slip channel, J. Comput. Math. 15 (1997) 335. [4] W. Bao, H. Han, High-order local artificial boundary conditions for problems in unbounded domains, Comput. Methods Appl. Mech. Engrg. 188 (2000) 455–471. [5] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [6] P.G. Ciarlet, the Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [7] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (1977) 629–651. [8] K. Feng, Asymptotic radiation conditions for reduced wave equations, J. Comput. Math. 2 (1984) 130–138.
H. Han, C. Zheng / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2011–2027
2027
[9] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992. [10] D. Givoli, J.B. Keller, Special finite elements for use with high-order boundary conditions, Comput. Methods Appl. Mech. Engrg. 119 (1994) 199–213. [11] D. Givoli, I. Patlashenko, An optimal high-order non-reflecting finite element scheme for wave scattering problems, Comput. Methods Appl. Mech. Engrg., to appear. [12] D. Givoli, I. Patlashenko, J. Keller, High-order boundary conditions and finite elements for infinite domains, Comput. Methods Appl. Mech. Engrg. 143 (1997) 13–39. [13] D. Givoli, I. Patlashenko, J. Keller, Discrete Dirichlet-to-Neumann maps for unbounded domains, Comput. Methods Appl. Mech. Engrg. 164 (1998) 173–185. [14] C.I. Goldstein, A finite element method for solving Helmholtz type equations in waveguides and other unbounded domains, Math. Comp. 39 (1982) 309–324. [15] T.M. Hagstrom, H.B. Keller, Exact boundary conditions at artificial boundary for partial differential equations in cylinders, SIAM J. Math. Anal. 17 (1986) 322–341. [16] T.M. Hagstrom, H.B. Keller, Asymptotic boundary conditions and numerical methods for nonlinear elliptic problems on unbounded domains, Math. Comp. 48 (1987) 449–470. [17] L. Halpern, M. Schatzman, Artificial boundary conditions for incompressible viscous flows, SIAM J. Math. Anal. 20 (1989) 308– 353. [18] H. Han, W. Bao, An artificial boundary condition for the incompressible viscous flows in a no-slip channel, J. Comput. Math. 13 (1995) 51–63. [19] H. Han, W. Bao, The artificial boundary conditions for incompressible viscous flows in channels, Lecture Notes Numer. Appl. Anal. 14 (1995) 47. [20] H. Han, W. Bao, An artificial boundary condition for two-dimensional incompressible viscous flows using the method of lines, Int. J. Numer. Methods Fluids 22 (1996) 483–493. [21] H. Han, J. Lu, W. Bao, A discrete artificial boundary condition for steady incompressible viscous flows in a no-slip channel using a fast iterative method, J. Comput. Phys. 114 (1994) 201–208. [22] H. Han, W. Bao, T. Wang, Numerical simulation for the problem of infinite elastic foundation, Comput. Methods Appl. Mech. Engrg. 147 (1997) 369–385. [23] H. Han, X. Wu, Approximation of infinite boundary condition and its application to finite element method, J. Comput. Math. 3 (1985) 179–192. [24] H. Han, X. Wu, The approximation of exact boundary condition at an artificial boundary for linear elastic equation and its application, Math. Comp. 59 (1992) 21–27. [25] J.L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Springer, Berlin, 1971. [26] D. Yu, Approximation of boundary conditions at infinity for a harmonic equation, J. Comput. Math. 3 (1985) 220–227.