The Riemann complex boundary element method for
the solutions of two-dimensional Elliptic equations
D.L. Young
*,T.J. Chang
1,T.I. Eldho
2Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, Taipei 10617, Taiwan Received 4 July 2000; received in revised form 2 January 2002; accepted 1 March 2002
Abstract
In this paper,a new boundary integral equation model Riemann complex boundary element method (RCBEM),is proposed based on the boundary element method (BEM) and the theory of Vekua and its modification as well as complex Riemann function as the fundamental solution. The RCBEM method is used to solve the linear,second order,elliptical partial differential equations in the fluid flow problems. In comparison to the generally used BEM,for RCBEM,there are two distinct differences. First one is that, RCBEM applies complex Riemann function as the fundamental solution of the adjoint operator while in direct BEM,on the other hand Green function is used. The second one is that the governing equations should be transformed into complex domain because there exist two characteristics in complex plane for elliptic systems,while in the direct BEM is not,since the Green function is adopted instead. The singular problem occurring in direct BEM can be avoided in RCBEM,especially for regular domain problems. The efficiency and accuracy of the RCBEM depends on the complex variable integration. To verify the feasi-bility and accuracy of RCBEM,the model is applied to different case studies of potential flows,Helmholtz equation problem and advection–diffusion problem and results are compared with analytical solutions and other numerical models. The results are satisfactory and prove the applicability of RCBEM for various two-dimensional elliptic equation problems.
Ó 2002 Elsevier Science Inc. All rights reserved.
Keywords: Riemann complex functions; Boundary element method; Fluid flow problems
www.elsevier.com/locate/apm
*
Corresponding author. Tel./fax: +886-2-23626114. E-mail address:dlyoung@hy.ntu.edu.tw(D.L. Young).
1
Present address: Department of Agricultural Engineering,National Taiwan University.
2
Present address: Department of Civil Engineering,IIT Bombay,India. 0307-904X/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 0 2 ) 0 0 0 4 8 - 3
1. Introduction
Boundary element method (BEM) has been established as a powerful numerical tool in the solution of various fluid flow problems [1]. In BEM,the computational domain becomes the enclosing boundary and the effective dimensions of the problem considered will be reduced by one. Hence it is much easier in discretization and data preparation for the problem considered. These advantages make BEM more suitable for the solution of various fluids flow problems.
BEM are usually derived from the Green’s theorem with an appropriate free-space Green’s function [2] using real variables. In recent times,a boundary element approach using complex variables for boundary integration,known as complex variable boundary element method (CVBEM) has been introduced [3]. The CVBEM is a generalization of the Cauchy integral for-mula into a boundary integral equation method,followed by forfor-mulation into a workable computer algorithm for effective mathematical simulation. This generalization allows an imme-diate and valuable transfer of the modeling technique and makes the process simpler and more efficient than using the real variables. But this limits its application to two-dimensional harmonic (Laplace) problems [3].
Generally an integral equation is solved by a numerical model that assumes the boundary of the problem domain is discretized into piecewise-polynomial curves,and the known and the unknown boundary values are approximated as piecewise-continuous functions along the boundary. Unlike the Green’s function formulations,the complex variable method does not depend on the shape of the contour between nodes [3,4]. Most of the complex variable boundary element methods adopt piecewise-linear representations of the complex functions that result in a second order accurate integration and normally give second-order accurate solutions for the boundary element solutions as well [4]. In complex BEM,the simplicity and elegance of complex analysis carries over to the computations as well. While the Green’s function formulation and the complex method are not directly compared here,studies by Dold and Peregrine [5] and Hromadka and Lai [3] indicate that the latter method is clearly superior to others in many fluid flow problems.
In the present study,a new complex boundary element method called Riemann Complex Boundary Element Method (RCBEM) is proposed. In RCBEM,the solutions of two-dimensional fluid flow problems are developed based on the theory of Vekua [6] and it’s modification and, thereby Riemann function [7] is taken as the fundamental solution to solve the linear second order elliptic equations.
In comparison with the generally used direct BEM,there are two major differences in RCBEM. In direct BEM,Greens functions are employed as solutions of adjoint operator while in RCBEM complex Riemann functions are used as the fundamental solutions of the adjoint operators. Secondly in RCBEM,the governing equations must be transformed into a complex domain be-cause there exist two characteristics in complex plane for elliptic systems (it should be noted that hyperbolic systems remain in real plane),while the direct BEM is solved in real planes,since the Green function is adopted instead. The main difference of RCBEM with the CVBEM [3] is,while in RCBEM,the Riemann functions are used as the fundamental solutions,in CVBEM,Cauchy’s functions are used as the fundamental solutions which will restrict the applications of more en-gineering problems beyond the harmonic functions.
In RCBEM,since the generic meaning of Riemann function is a characteristic boundary value problem,and the characteristic curves are closely associated with the propagation of certain types
of singularities,we can infer that the Riemann function used is a regular solution,meaning that the singularities are far behind infinite. Some of the important features offered by RCBEM are: (1) the Riemann functions used are regular solutions and hence they satisfy the governing equations throughout the region enclosed by the problem boundary,the approximation is made only at the boundary; (2) the integration of the boundary integrals along each boundary element are carried out exactly; (3) mathematical means can be devised to evaluate approximation errors; and (4) substantial modeling simplification are possible resulting from the complex variable application and the boundary element approach. However,it should be noted that the present theory of RCBEM is for two-dimensional problems only.
In this paper,the RCBEM is applied for the solutions of various fluid flow problems governed by Laplace equation,Helmholtz equation and steady state advection–diffusion equations. The feasibility and accuracy of the RCBEM has been shown by solving a variety of problems governed by the above mentioned equations. RCBEM solutions are verified by comparing with available exact solutions and direct BEM solutions. Good agreements are observed in all the cases.
2. Complex Riemann boundary integral equations 2.1. Riemann function in complex variables
Initially,let us take account of some results of the theory of solutions of linear second order elliptic differential equations essentially based on the complex Riemann function and the Volterra type integral equation described by Vekua [6]. Consider a second order linear elliptic differential equation in two independent variables x and y,
L½u ¼o 2u ox2þ o2u oy2þ aðx; yÞ ou oxþ bðx; yÞ ou oyþ cðx; yÞu ¼ f ðx; yÞ ð2:1Þ
where the coefficients a; b; c and f are functions of the variables x and y,are analytic and con-tinuous in a domain D. Let the coordinate transformation of the complex variables be represented as z¼ x þ iy, zz ¼ x iy and the differential operators as,
o oz¼ 1 2 o ox i o oy ; o ozz¼ 1 2 o ox þ i o oy
Now using the new complex variables,Eq. (2.1) can be transformed into, M½u ¼ o 2u ozozzþ Aðz; zzÞ ou ozþ Bðz; zzÞ ou ozzþ Cðz; zzÞu ¼ F ðz; zzÞ ð2:2Þ
which is a complex form of hyperbolic partial differential equation. A, B, C and F are defined as,
Aðz; zzÞ ¼1 4 a zþ zz 2 ; z zz 2i ! " þ ib zþ zz 2 ; z zz 2i !# ; Cðz; zzÞ ¼1 4c zþ zz 2 ; z zz 2i !
Bðz; zzÞ ¼1 4 a zþ zz 2 ; z zz 2i ! " ib zþ zz 2 ; z zz 2i !# ; Fðz; zzÞ ¼1 4f zþ zz 2 ; z zz 2i !
Eq. (2.2) is the complex form of (2.1),and the coefficients A, B, C and F are holomorphic with respect to variables z and zz in a domain z2 D and zz 2 D1 where D and D1 are simply connected domains.
Now the Riemann function can be expressed in two alternative forms [6], Rðt; zz; t; sÞ ¼ exp Z zz s Aðt; gÞ dg " # on z¼ t ð2:3Þ Rðz; s; t; sÞ ¼ exp Z z t Bðn; sÞ dn on zz¼ s ð2:4Þ
where t and s are fixed parameters. The Riemann function R satisfies the following normalized condition,
Rðz; s; t; sÞ ¼ 1 on z¼ t; zz ¼ s ð2:5Þ
By taking the adjoint of Eq. (2.2),using the Riemann functions and integrating with respect z, zz, we have the Volterra type integral equation of second kind [6],
Rðz; zz; t; sÞ Z zz s Aðz; gÞRðz; g; t; sÞ dg Z z t Bðn; zzÞRðn; zz; t; sÞ dn þ Z z t Z zz s Cðn; gÞRðn; g; t; sÞ dn dg ¼ 1 ð2:6Þ
As described in Vekua [6],using the adjoint property of Eq. (2.2) and using the Riemann func-tions,exchanging the pairs of ðz; zzÞ and ðt; sÞ and integrating with respect to t in the interval ðz0; zÞ; s in ðzz0; zzÞ,one can obtain
uðz; zzÞ ¼ uðz0; zz0ÞRðz; zz; z0; zz0Þ þ Z z z0 U1ðtÞRðz; zz; t; zz0Þ dt þ Z zz zz0 U2ðsÞRðz; zz; z0;sÞ ds þ Z z z0 Z zz zz0 Fðt; sÞRðz; zz; t; sÞ ds dt ð2:7Þ where U1ðzÞ ¼ ouðz; zz0Þ oz þ Bðz; zz0Þuðz; zz0Þ; U2ðzÞ ¼ ouðz0; zzÞ
ozz þ Aðz0; zzÞuðz0; zzÞ.
It should be noted that the coefficients of Eq. (2.1) are real functions and u is also a real function and hence after taking integration by parts,the solution of (2.1) can be written as,
uðx; yÞ ¼ Re H0ðzÞ/ðzÞ " þ Z z z0 Hðz; tÞ/ðtÞ dt þ Z z z0 Z zz zz0 Rðz; zz; t; sÞF ðt; sÞds dt # ð2:8Þ
where Re means real part and
Hðz; tÞ ¼ o
otRðz; zz; t; z0Þ þ Bðz; zz0ÞRðz; zz; t; z0Þ ð2:10Þ /ðzÞ ¼ 2uðz; zzÞ uðz0; zz0ÞRðz; zz0; z0; zz0Þ ð2:11Þ /ðzÞ is an arbitrary holomorphic function in D determined by the boundary condition. Assuming without loss of generality that z0 ¼ 0 lying inside the region D and also Að0; zzÞ ¼ Bðz; 0Þ ¼ 0,we get,
Hðz; tÞ ¼ o
otRðz; zz; t; 0Þ ð2:12Þ
/ðzÞ ¼ 2uðz; 0Þ uð0; 0ÞRðz; 0; 0; 0Þ ð2:13Þ
The detailed procedure of the derivation of Eq. (2.8) is given in Vekua [6],which is not repeated here.
2.2. Riemann functions for the elliptic differential equations
Before the application of boundary element procedure for the concerned elliptic differential equations using RCBEM,we have to find the fundamental solutions (Riemann functions in RCBEM). The fundamental solution is derived from the general Eq. (2.6). In this section,we consider some important elliptic equations in different forms,such as Laplace equation,modified Helmholtz equation and Helmholtz equation. Equations like steady-state advection–diffusion equation are solved after converting it into the modified Helmholtz equation.
2.2.1. Laplace equation
Considering the Laplace equation,
r2u¼ 0 ð2:14Þ
Comparing with the Eq. (2.2),the coefficients Aðz; gÞ ¼ Bðn; zzÞ ¼ Cðn; gÞ ¼ 0. From Eq. (2.6),the fundamental solution (Riemann function) for the Laplace equation is obviously,
Rðz; zz; t; sÞ ¼ 1 ð2:15aÞ
Hence the Eqs. (2.9) and (2.12) for potential flow problems can be written as,
H0ðzÞ ¼ 1 ð2:15bÞ
Hðz; tÞ ¼ o
otRðz; zz; t; oÞ ¼ 0 ð2:15cÞ
2.2.2. Modified Helmholtz equation
Consider the modified Helmholtz equation,
r2u k2u¼ 0 ð2:16Þ
where k is assumed a constant. In comparison with Eq. (2.2),the coefficients Aðz; gÞ ¼ Bðn; zzÞ ¼ 0, and Cðn; gÞ ¼ k2=4. Now Eq. (2.6) can be expressed as,
Rðz; zz; t; sÞ 1 4k 2 Z z t Z zz s Rðn; g; t; sÞ dn dg ¼ 1 ð2:17Þ
By the methods of successive approximations [8],the solution is, Rðz; zz; t; sÞ ¼ 1 þ Z z t Z zz s C0ðz; zz; n; g; t; sÞ dg dn ð2:18Þ where C0ðz; zz; n; g; t; sÞ ¼P1m¼1N0ðmÞðz; zz; n; g; t; sÞ N0ðmÞðz; zz; n; g; t; sÞ ¼ Z z n Z zz g N0ðm1Þðs; r; n; g; t; sÞN0dr ds; N0ð1Þ¼ N0 ¼ Cðn; gÞ since N0ð1Þ¼ Cðn; gÞ ¼ k2=4; Nð2Þ 0 ¼ Z z n Z zz g 1 4k 2 1 4k 2 dr ds¼k 4 16ðz nÞðzz gÞ N0ð3Þ¼ Z z n Z zz g 1 4k 2 k4 16ðz nÞðzz gÞ dr ds¼k 6 64ðz nÞ 2 ðzz gÞ2. . . Now we have, Rðz; zz; t; sÞ ¼ 1 þ Z z t Z zz s C0ðz; zz; n; g; t; sÞ dg dn ¼ 1 þ Z z t Z zz s X1 m¼1 N0ðmÞðz; zz; n; g; t; sÞ " # dg dn ¼ 1 þ Z z t Z zz s 1 4k 2 þk 4 16ðz nÞðzz gÞ þ k6 64ðz nÞ 2 ðzz gÞ2þ dg dn ¼X 1 k¼1 ½ðz tÞðzz sÞ=22k=k!Cðk þ 1Þ ¼ I0 k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz tÞðzz sÞ p ð2:19aÞ where I0is the modified Bessel function of the first kind of order zero which also shows the regular behavior in the defined domain.
Hence the Eqs. (2.9) and (2.12) can be written as, H0ðzÞ ¼ I0ðkrÞ; r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2þ y2 p ð2:19bÞ Hðz; tÞ ¼1 2k ffiffi zz p I1 k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi zzðz tÞ p ffiffiffiffiffiffiffiffiffiffi z t p ð2:19cÞ
where Ij is the modified Bessel function of the first kind of order one which is a regular function for the defined domain.
2.2.3. Helmholtz equation
For the Helmholtz equation (k is constant), r2uþ k2
In comparison with Eq. (2.2),the coefficients Aðz; gÞ ¼ Bðn; zzÞ ¼ 0,and Cðn; gÞ ¼ k2=4. Using the successive approximation as mentioned above,we obtain
Rðz; zz; t; sÞ ¼ J0 k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz tÞðzz s p Þ ð2:21aÞ where J0 is the Bessel function of the first kind of order zero,which is also a regular function for the defined domain.
Hence the Eqs. (2.9) and (2.12) can be written as, H0ðzÞ ¼ J0ðkrÞ; r¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2þ y2 p ð2:21bÞ Hðz; tÞ ¼ o otRðz; zz; t; 0Þ ¼ 1 2k ffiffi zz p J1 k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi zzðz tÞ p ðpffiffiffiffiffiffiffiffiffiffiz tÞ ð2:21cÞ
where J1 is the Bessel function of the first kind of order one,which shows the regular behavior in the defined domain.
2.3. Boundary conditions
In this section,the application of two types of boundary conditions,Dirichlet boundary con-ditions and mixed boundary concon-ditions,are illustrated with reference to the RCBEM model development.
2.3.1. Dirichlet boundary conditions
Let D be a simply connected domain bounded by a contour C. Considering the function uðx; yÞ satisfying the real boundary condition,
uðx; yÞ ¼ gðx; yÞ ¼ gðsÞ s2 C ð2:22Þ
Assuming that gðsÞ satisfies the H€oolder condition on boundary C and the boundary value /ðtÞ of Eq. (2.13) also satisfies the H€oolder condition which can be expressed in the form of double layer (or the Cauchy kernel integral),
/ðzÞ ¼ 1 2pi Z C /ðtÞ t zdt ð2:23Þ
for harmonic function [6] and /ðzÞ ¼ 1 2pi Z C lðtÞ t zdt ð2:24Þ
for the other functions [6] where lðtÞ is a real density function of the dipole.
Substituting Eq. (2.24) into Eq. (2.8),the Fredholm integral equation of the first kind can be obtained as,
uðx; yÞ ¼ Z
C
where Mðz; tÞ ¼ Re H0ðzÞdt=ds 2piðt zÞ þdt=ds 2pi Z z 0 Hðz; t1Þ t t1 dt1 ð2:26Þ Let an interior point z approaches to an arbitrary point t0,and satisfies the boundary condition Eq. (2.22),the Fredholm type integral equation of second kind can be obtained,
1 aa 2p ! lðt0Þ þ Z C Mðt0; tÞlðtÞ ds ¼ gðt0Þ ð2:27Þ
for the unknown function l, (aa is contour angle,0 6 aa 62pÞ [6],where, Mðt0; tÞ ¼ Re H0ðt0Þdt=ds 2piðt t0Þ þdt=ds 2pi Z t0 0 Hðt0; t1Þ t t1 dt1 ð2:28Þ Eqs. (2.25)–(2.28) are convenient to evaluate in numerical methods since the integral equations are all regular.
2.3.2. Mixed boundary conditions
Considering the following more general mixed type boundary condition, pðsÞou
oxþ qðsÞ ou
oyþ rðsÞu ¼ gðsÞ s 2 C ð2:29Þ
and the fact that o ox¼ o ozþ o ozz; o oy ¼ i o oz o ozz ; Eq. (2.29) can be written as
AðsÞou ozþ AðsÞ
ou
ozzþ rðsÞu ¼ gðsÞ ð2:30Þ
where AðsÞ ¼ pðsÞ þ iqðsÞ, AðsÞ ¼ pðsÞ iqðsÞ. Substituting Eq. (2.30) into Eq. (2.8) yields,
Re Lðt0Þ/0ðt0Þ þ Mðt0Þ/ðt0Þ þ Z t0 0 Nðt0; t1Þ/ðt1Þ dt ¼ gðt0Þ ð2:31Þ
where Lðt0Þ ¼ Aðt0ÞH0ðt0Þ and
Mðt0Þ ¼ Aðt0Þ oH0ðt0Þ ot0 þ Aðt0ÞH ðt0; tÞ þ Aðt0Þ oH0ðt0Þ ott0 þ rðt0ÞH0ðt0Þ ð2:32Þ Nðt0; t1Þ ¼ Aðt0Þ oH0ðt0; t1Þ ot0 þ Aðt0Þ oHðt0; t1Þ ott0 þ rðt0ÞH ðt0; t1Þ ð2:33Þ
Assuming that the boundary values /ðzÞ and /0ðzÞ satisfy the H€oolder condition on boundary C, we can use an integral representation of the following,
/ðzÞ ¼ 1 2pi
Z C
lðtÞ lnðt zÞ dt z2 D ð2:34Þ
With the help of Eqs. (2.34),(2.31) can be written as, Re ( 1 aa 2p ! Lðt0Þlðt0Þ Z C Lðt0Þ t0 t E0ðt0; tÞ lðtÞ dt ) ¼ gðt0Þ ð2:35Þ where, Eðt0; tÞ ¼ 1 2pi Z t0 0 Nðt0; tÞ lnðt t1Þ dt1þ Mðt0Þ lnðt t0Þ ð2:36Þ Eqs. (2.34)–(2.36) are convenient to evaluate in numerical methods since the integral equations are all regular. The above equations can be used in the evaluation of mixed type and Neuman type boundary conditions in the RCBEM procedure.
3. Boundary integral procedure
Using the results described in Section 2,RCBEM is developed to solve the singular integral equation. The first step in the RCBEM procedure is that the boundary of the domain of the problem considered is divided into piecewise smooth elements as shown in Fig. 1. On any one of the elements considered,the boundary conditions are of the same kind,say Dirichlet,Neumann or mixed etc. Assuming t0 as a base point (as in Fig. 1),we have
t t0 ¼ reia ð3:1Þ
dt¼ eihds ð3:2Þ
Then the kernel of Eq. (2.27) or (2.28) can be expressed as, Mðt0; tÞ ¼ 1 2pIm H0ðt0ÞeiðhaÞ r þ eih Z t0 0 Hðt0; t1Þ t t1 dt1 ¼ 1 2p H0ðt0Þ sinðh aÞ r þ 1 2pIm e ih Z t0 0 Hðt0; t1Þ t t1 dt1 ð3:3Þ The first term of Eq. (3.3) is regular as t approaches to t0 [8].
Our aim is to finally devise a method that will convert Eq. (2.27) into a series of algebraic equations. For convenience,considering a linear element,the unknowns of Eq. (2.27) can be represented as follows [9], lðtÞ ¼X 2 i¼1 Nili¼ ðljþ1 ljÞn þ njþ1lj njljþ1 ðnjþ1 njÞ nj6n 6 njþ1 ð3:4Þ
in which n is the distance along the element (see Fig. 2), Ni is the shape function of the linear element and, ri ¼ ðn2þ g2iÞ 1=2 ð3:5Þ ai ¼ tan1 nsin hþ gicos h ncos h gisin h ð3:6Þ are used in the integration.
The last term of Eq. (3.3) is the imaginary part of complex variable integration,which can be evaluated using the simple trapezoidal rule,Simpson 3/8 rule or by Gaussian quadrature. Here the implementation using Gaussian quadrature is explained briefly.
In the Gaussian quadrature method,the complex variable integration is done by separating the real and imaginary parts,and integrating each separately. For example,consider the Helmholtz equation, Hðt0; t1Þ ¼ k 2 ffiffiffiffi tt0 p J1 k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tt0ðt0 t1Þ p ffiffiffiffiffiffiffiffiffiffiffiffiffi t0 t1 p ð3:7Þ
Substituting Eq. (3.7) into Eq. (2.28),the complex integration can be expressed as, Z t0 0 J1ðk ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tt0ðt0 t1Þ p Þ ðtj t1Þ ffiffiffiffiffiffiffiffiffiffiffiffiffit0 t1 p dt1 0 6 t16t0; t0 ¼ x0þ iy0 ð3:8Þ
Since Eq. (3.8) is independent of path,a straight-line integration can be chosen between 0 and t0. Let t1ðsÞ xðsÞ þ iyðsÞ ¼ s þ i y0 x0 s 0 6 s 6 x0 ð3:9Þ dt1 ¼ 1 þ iy0 x0 ds ð3:10Þ
Substituting Eqs. (3.9) and (3.10) into Eq. (3.8),the integral can be separated into real part and imaginary part and the Gaussian quadrature can be applied for numerical integration.
The last term of Eq. (3.3) can be written as follows after the numerical discretization for each li 1 2pIm e ih Z t0 0 Hðt0; t1ÞNj ti t1 dt1 lj ¼ Bijlj; i¼ 1; 2; . . . N ð3:11Þ
The integral on the first term of Eq. (3.3) can be integrated over the element between pj and pjþ1 with respect to the base point pi,by omitting the constant of H0ðt0Þ.
Ie ¼ Z njþ1 nj sinðhj aiÞ ri lðnÞ dn ¼ ½Ke lj ljþi ð3:12Þ where jKej ¼ jðKeÞ i;j;ðK eÞ i;jþ1j ¼ jnjþ1I11 I12; I12 njI11j ð3:13Þ in which I11¼ ðsin hiI1 cos hjI2Þ 1 2p; I12¼ ðsin hjI3 cos hjI4Þ 1 2p ð3:14Þ I1 ¼ 1 ðnjþ1 njÞ Z njþ1 nj cos ai ri dn¼ 1 ðnjþ1 njÞ 1 2 cos hjln n2jþ1þ g2 i n2j þ g2 i ! " sin hj tan1 njþ1 gi tan1 nj gi # ð3:15Þ I2 ¼ 1 ðnjþ1 njÞ Z njþ1 nj sin ai ri dn¼ 1 ðnjþ1 njÞ 1 2sin hjln n2jþ1þ g2 i n2jþ g2 i ! " cos hj tan1 njþ1 gi tan1 nj gi # ð3:16Þ
I3¼ 1 ðnjþ1 njÞ Z njþ1 nj cos ai ri ndn ¼ 1 ðnjþ1 njÞ cos hj ðnjþ1 ( njÞ gi tan1 njþ1 gi tan1nj gi gj 2 sin hjln n2jþ1þ g2 i n2jþ g2 i !) ð3:17Þ I4¼ 1 ðnjþ1 njÞ Z njþ1 nj sin ai ri ndn ¼ 1 ðnjþ1 njÞ sin hj ðnjþ1 ( njÞ gi tan1 njþ1 gi tan1nj gi þgi 2 cos hjln n2jþ1þ g2 i n2jþ g2 i !) ð3:18Þ Finally,the algebraic Eq. (2.27) can be expressed as,
XN j¼1 ðAijþ BijÞlj¼ gi i¼ 1; 2; . . . N ð3:19Þ where Aij¼ ðKeÞij þ 1 a a 2p ! dij dij¼ 1 i¼ j 0 i6¼ j ð3:20Þ Eq. (3.20) is used only in the case of interior problems. For the exterior problems,Eq. (3.20) is written as, Aij¼ ðKeÞij a a 2pdij dij ¼ 1 i¼ j 0 i6¼ j ð3:21Þ It should be noted that for exterior problem,the normal vector is in opposite direction.
4. Numerical applications
To test the feasibility and accuracy of the RCBEM model,here three numerical applications are presented. Initially,we consider the potential flow problems in which interior and exterior problems are considered. Secondly,a fluid flow problem governed by the Helmholtz equation is considered and finally a steady-state advection–diffusion problem is solved by transforming the governing equation into the modified Helmholtz equation. The RCBEM results in all the cases are compared with analytical and other numerical solutions.
4.1. Potential flow problems
The governing equation for the potential flow problems is the Laplace equation. Here the potential flows are considered with interior problems and exterior problems. In the case of
po-tential flow problems,we consider the solution of two problems using the RCBEM. Initially an interior problem of groundwater flow is considered and then an exterior problem of source outside a circular cylinder is considered.
4.1.1. Groundwater flow problem
Here we consider a two-dimensional groundwater flow problem in a homogeneous,isotropic porous media with a pumping well. The problem is a bounded square region,0 6 x 6 1 and 0 6 y 6 1 with zero potential on all the sides as shown in Fig. 3. There is a pumping well at the middle of the domain with unit discharge. In the RCBEM model,the boundary of the domain is discretized with 40 linear elements. Fig. 4 shows the potential variation with respect to the x-axis at y¼ 0:5 from the well position. The results are compared with exact solution and a BEM model [2]. Good agreement is observed between the results. A sensitivity study with different mesh discretization showed that RCBEM yields comparable results with other models even with 8 linear boundary elements. This case study shows the feasibility of the RCBEM model for the interior type problems.
4.1.2. Exterior flow problem
Here we consider the two-dimensional flow field past a circular cylinder of radius r ¼ 1 due to a source at x¼ 2 with strength m ¼ 2p. The aim is to find the stream function distribution over the circular cylinder due to the source effects. The problem domain with discretization is shown in Fig. 5. The boundary conditions of the stream function along the cylinder are assumed to be zero. Since the problem is exterior in nature,a sensitivity study of the mesh showed that a finer dis-cretization is necessary [2]. In the RCBEM model,the total boundary of the domain is discretized with 160 linear elements. The problem can be simplified by considering the symmetry and dis-cretizing half of the domain.
In this problem,the influence of the point source can be directly taken into account or it can be separated into the perturbed part and the potential flow part and the solution can be combined [2]. In the present analysis,as the problem considered is linear in nature,it is convenient to separate the stream function (w) into two parts, w¼ wUþ wP,where wUdefines the stream function due to potential flow and wP is a perturbed stream function owing to the influence of the source. Fig. 6 shows the stream function variation around the circular cylinder. The results are compared with exact solution. Good agreement is observed between the results. This case study shows the fea-sibility of the RCBEM model for the exterior type problems.
Fig. 4. Comparison of potential for groundwater flow problem.
4.2. Helmholtz equation problem
The governing equation of the reduced wave problem is the Helmholtz equation. To illustrate the application of the Helmholtz equation problem,here we consider an interior problem in which the potential is propagated from a point source of unit strength. For the demonstration pur-pose,the value of k is assumed as 1. The problem domain with discretization is shown in Fig. 7. The point source is located at the center of the domain. All over the boundary of the domain,a Dirichlet boundary condition uð1; hÞ ¼ 1 is assumed. The boundary of the domain is discretized with 120 elements. The exact solution for the problem considered is given by,
u¼J0ðkrÞ J0ðkaÞ
ð4:1Þ Fig. 8 shows the propagation of potential for a¼ 1 along the axis using the RCBEM. The results are compared with the exact solution and a BEM model [2]. The results are almost identical. A sensitivity study with different mesh discretization showed that RCBEM yields
Fig. 6. Streamlines for the exterior potential flow problem: comparison of RCBEM and exact solution.
comparable results with other models even with 40 linear boundary elements. This case study shows the feasibility of the RCBEM model for the Helmholtz equation type problems.
4.3. Advection–diffusion problem
Here the steady-state advection–diffusion problem is solved,by converting the governing equation into the modified Helmholtz equation. The governing equation of steady-state advec-tion–diffusion problem is,
r ð½kruÞ ðV rÞu ¼ f ð4:2Þ
where k is the diffusivity coefficient, V is the convective velocity and f is the source or sink term. The governing equation is made dimensionless using the following variables,
xx ¼xx L; r ¼ Lr; ½k ¼½k K ; V ¼V U; u ¼ u u0 ð4:3Þ where L is the characteristic length, K is the characteristic diffusivity and U is the characteristic velocity. Assuming f ¼ 0,the non-dimensional form of Eq. (4.2) can be written as,
r ð½kruÞ PeðV rÞu ¼ 0 ð4:4Þ
where Pe¼ UL=K is the Peclet number. Using the following transformation [10], u¼ egu; rg¼1
2Pe½k 1
V ð4:5Þ
into Eq. (4.4),we can obtain the transformed modified Helmholtz equation as the following equation after omitting the superscript for brevity,
r ð½kruÞ k2u¼ 0 ð4:6Þ where k2¼1 4Pe 2V½k1 V þ1 2Per V ð4:7Þ
There are two limitations for this transformation viz. k2 have to be positive and the flow field should be irrotational. Here we consider an advection–diffusion problem in two-dimensional domain. The problem is a bounded square domain. A boundary condition with sinðpyÞ concen-tration profile is imposed on the boundaryð0; yÞ and zero concentration is assumed on all other boundaries. Fig. 9 depicts the studied domain with associated boundary conditions. The RCBEM simulates the flow field with constant Vx¼ 1, Vy ¼ 0, ½k ¼ ½I,and different Peclet numbers from diffusion (low Pe) to advection (high Pe) dominated flows. The boundary of the domain is dis-cretized into 40 linear elements. The computations are carried out for Peclet numbers of 0,1,10, 40 and 80. For this problem,an exact solution can be derived as,
u¼ sin py ea eb½e aþbx eaxþb ð4:8Þ where a¼Peþ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4p2þ Pe2 p 2 ; b¼ Pepffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi4p2þ Pe2 2 ð4:9Þ
Using RCBEM,for all the Peclet numbers,very stable and comparable results with analytical solutions are obtained. Fig. 10 shows the concentration along ðx; 0:5Þ compared with exact so-lution for various Peclet numbers. Good agreement is observed between the soso-lutions and the results show that RCBEM is stable for higher Peclet number problems,which are more difficult and challengable to other numerical methods. This case study shows the effectiveness of the RCBEM model for the two-dimensional steady-state advection–diffusion problems.
5. Concluding remarks
A new boundary integral equation model,RCBEM,is proposed based on the (BEM) and the theory of Vekua and its modification as well as complex Riemann function as the fundamental solution. The RCBEM method is used to solve the linear,second order,elliptical partial differ-ential equations in the fluid flow problems. Compared to the generally used BEM,RCBEM applies complex Riemann function as the solution of the adjoint operator and the governing equations should be transformed into complex domain. The singular problem occurring in direct BEM can be avoided in RCBEM,especially for regular domain problems. The efficiency and accuracy of the RCBEM depends on the complex variable integration and the presented meth-odology applies only to two-dimensional problems. The simulations of various numerical ex-amples presented demonstrated the accuracy and feasibility of RCBEM model.
Acknowledgements
The work reported in this paper was supported by the National Science Council,Taiwan. It is greatly appreciated.
References
[1] H. Power,L.C. Wrobel,Boundary Integral Methods in Fluid Mechanics,Computational Mechanics Publications, Southampton,1995.
[2] C.A. Brebbia,J.C.F. Telles,L.C. Wrobel,Boundary Element Techniques––Theory and Applications in Engi-neering,Springer,Berlin,1984.
[3] T.V. Hromadka II,C. Lai,The Complex Variable Boundary Element Method in Engineering Analysis,Spinger, New York,1986.
[4] W.W. Schultz,S.W. Hong,Solution of potential problems using an over determined complex boundary integral method,J. Comp. Phys. 84 (1988) 414–440.
[5] J.W. Dold,D.H. Peregrine,Steep unsteady water waves: an efficient computational scheme,in: Proceedings of the 19th International Conference On Coastal Engineering,Houston,TX,USA,1984.
[6] I.N. Vekua,New Methods for Solving Elliptic Equations,North-Holland Publishing Company,Amsterdam,1967. [7] R. Courant,D. Hilbert,Methods of Mathematical Physics,vol. 1,Interscience,New York,1953.
[8] F.G. Tricomi,Integral Equations,Interscience,New York,1957.
[9] V.V. Ivanov,The Theory of Approximate Methods and their Application to the Numerical Solution of Singular Integral Equation,Nordhoff Publishing Co.,Leyden,1976.
[10] M. Ikeuchi,M. Sakakihara,Boundary elements in steady convective diffusion problems,J. Comput. Appl. Math. 12/13 (1985) 381–389.
