Boundary element analysis of the three-dimensional convective diffusion equation

Fusion Engineering and Design 19 (1992) 315-319 North-Holland

315

Boundary element analysis of the three-dimensional convective diffusion equation Toshihisa Honma

Department of Electrical Engineering, Faculty of Engineering, Hokkaido University, Sapporo, Japan Accepted 26 June 1992 Handling Editor: K. Miya

The present paper shows boundary element solutions for three-dimensional convective diffusion equations with constant velocity field. In particular, mixed boundary element solutions in steady- and unsteady-state field are presented, and the usefulness of the method is discussed.

Introduction In recent years, computational fluid dynamics is increasingly developed, and is giving a number of valuable results which cannot be obtained from experimental data. However, it is still very difficult to numerically solve a convective diffusion equation when the convective term dominates [1,2]. In addition, it is also difficult to solve a three-dimensional problem, even by using powerful numerical methods such as the finite difference method and the finite element method, though we fully recognize that three-dimensional analysis is a most significant problem in engineering analysis. In the present paper, the author shows boundary element solutions for three-dimensional convective diffusion equations with constant velocity field [3-5]. In particular, mixed boundary element solutions in steadyand unsteady-state fields are presented, and the usefulness of the method is discussed.

phenomena in classical physics with which we are familiar in various fields of science and engineering. Now we shall consider a bounded domain O enclosed by the surface F. It is supposed that the conservative law for some physical quantity ~b holds in this region:

a6 --+v.qfp,

(1)

at

in which p stands for the source/sink term, and the vector q denotes the total mass flux. In this configuration, we assumed that the flux q is constructed from the diffusive flux qd -----KV~b, based on Fick's first law with the scalar diffusivity K, and the convective flux ¢c = v~b, v being the convection velocity. Accordingly, from eq. (1), we may describe the convective diffusion equation with respect to ~b as

a6

--

at

+ v.(-Kv~)

+ v . v~ =p,

(2)

2. Convective diffusion equation

where we impose the continuity equation for incompressible, viscous flow: V- v = 0.

As is well-known, the convective diffusion equation is based upon the conservation law of mass in classical physics and so mathematically formulated one of the most basic governing equations describing macroscopic

3. Steady-state convective diffusion equation 3.1. Boundary element formulation

Correspondence to: Dr. Toshihisa Honma, Dept. of Electrical

At the outset, we shall treat a steady-state field, which may be expressed using eq. (2) by

Engineering, Faculty of Engineering, Hokkaido University, Nishi-8, Kita-13, Kita-Ku, Sapporo 060, Japan.

V ' ( - x V ~ b ) + v " V~ =p.

0 9 2 0 - 3 7 9 6 / 9 2 / $ 0 5 . 0 0 © 1992 - E l s e v i e r Science P u b l i s h e r s B.V. A l l rights r e s e r v e d

(3)

316

T. Honma / The convective diffusion equation

When we may regard the velocity field as a constant in space, then the boundary integral equation equivalent to eq. (3), may be given [3,6] by

c(ri)6(ri) -frp*(r; r i ) ' n c ~ ( r ) = -fr~*(r;

(a,O,e) V

q)ad(r)'n d r

r i ) = exp[ /

Ivllr'l 2K

]

/(47rK I r ' l ) ,

(5) p * ( r ; ri) = - K V ~ F * ( r ; ri) - v q r * ( r ; ri) ,

ri).n dr.

17 ,..- ...."'"'""'"'" ................ ~ a ~ b p O ) "' J~- X (0,0,0) (a,0,0) Fig. 1. Three-dimensional convective diffusion domain to be analyzed.

(6)

and from the knowledge of the conservative law, c(ri) = fp*(r;

Y

(4)

where we introduce the fundamental solution of the adjoint field associated with eq. (3) in three dimensions [6]: v'r'+

(a,b,e)

(O,b,e)

dr

+fp(r)~*(r)~l"*(r; ri) dO,

~*(r;

Z

(7)

Here, the diffusivity K is assumed to be uniform over the whole domain, the vector r stands for the position vector of a field point e ~ u F, r i denotes the position vector of the source point where we want to know the field, and so r ' - - r - r i. Discretization of the boundary of the domain in planar patches such as triangles or quadrangles, approximation of the dependent variable and its normal diffusive flux within a patch by employing for instance the mixed element [7] and numerical boundary integration of eq. (4), may lead to the standard (or mixed) boundary element matrix equation [3-6].

3.2. Boundary element solution In order to demonstrate the boundary element solution, as example, we shall treat a cubic domain as shown in fig. 1, where it is assumed that the vector v = (v~, 0, 0), vx = constant, and the dimensions a = b = c = 1 for simplicity. Figure 2 shows typical boundary element solutions in the case of a boundary mesh size h =1~, along the center line (x, y, z ) = (x, 7, 1 7). The boundary element solution (small circle) is compared with the exact solution (solid line). It is observed that stable and accurate solutions are obtained even at high P6clet number P e ( = l vii~K), l being a typical length of the region. Defining the cell P~clet number as

P¢,c = Iv I h/K, the result is equal to P¢,c= 3.3 and 16.7, which exceeds the stability criterion Pe,c ~< 2 in the standard finite difference method and finite element method [1,2]. Similar results are also obtained by using a constant element [6] and a regular boundary element [2].

o.6r

v

00,4 0

Colcu'c ted ExQct

0.2

0

I 0

I 0.2

I

I 0.4

I

I 0.6

I

I 0.8

I 1.0

x/O

Fig. 2. Steady m i x e d b o u n d a r y element solution. A D i r i c h i e t b o u n d a r y c o n d i t i o n is imposed: (#(0, y, z)=sin(~-y/b).

sin(~z/c),

@(a, y, z) = @(x, O, z) = &(x, b, z ) = 4)(x, y, 0) =

~(x, y, c) = O.

Z Honma / The convective diffus~n equation Cell Peclet Number 20

40

60

80

100

120

140

Pe,c

160

317

[Dimensionless]

180

200

220

240

260

280

300

320

340

2 1 0 te

-1 -2

h21/6

h=1/8

-3 -4

-S -6 -7 -8

-9 -10

Fig. 3. Dependence of the relative error on the cell Peclet number.

Furthermore, we tried to analyze a convective diffusion problem for larger P6clet numbers. Figure 3 shows the dependence of the relative error on the cell P6clet number [9]. It is found that the relative error even less depends upon the cell P4clet number of P6clet number, as the P6clet number becomes larger. However, it is also shown that the boundary element solution is not

.......o''" ..r"

....-'; Pe= 1 0

.... " i

1.o r •""

T=O

. . . .

!

......

\ e~o.s

I.......i J/ ! 0.6

L.~'" 1 .o~. I"i :

0.4

.....

~ =

01 2S

h=I/4

~

!......I

"~

unconditionally stable. This may be due to numerical integration of a singular boundary integral. It is expected that the result may be improved by developing a more precise numerical evaluation of the boundary integrals. As for a variable velocity field, some boundary element formulations are presented. Details of these

......

.b'"

!

10 - 2 ~¥~e

0.2

o

0.2

0.4

x/a

0.6

( y/b=Z/c=l/2

0.8

1.0

)

Fig. 4. Unsteady mixed boundary element solution. Homogeneous Dirichlet boundary conditions and initial condition ~b0 = sin(~-x / a ). sin(~r y / b)- sin(~rz / c ) are imposed.

T. Honma / The convective diffusion equation

318

methods will not be given here but can be found elsewhere [10-2].

4. Unsteady-state convective diffusion equation

4.1. Boundary element formulation Next, we shall consider the unsteady-state convective diffusion equation (2). When the velocity is not a function of both time and space, then we may derive a boundary integral equation equivalent to eq. (2) [4,5,13]: :'r

/-

jo J/*(r, t r,. =_

:o7 qt*(r, t l r i , , r ) q d ( r , t ) . n d F

+ f ~'*(r, 01r~, ~)~o(r) +

n (r, ,)d. d,

(8)

t l r i, 7)

[

v'r'

Ivlet '

Ir'l 2

- (4~rKt,)a/2 exp[

2K

4K

4Kt'

Hit'] [ (47rKt,)d/2 exp

[

v'r' 2K

Ir'l ~2Kt'

(lO)

dO

~ * ( r , t l r i, r ) p ( r , t) dO,

n[t']

~ * ( r , t[ri,'r )

dt

where we utilize the time-dependent fundamental solution of the adjoint equation in d-dimensional space [13] ~*(r,

the convective and diffusive term in the singular kernel functions of the integrands. Hence, the integration must be numerically performed, but we have to expend much cpu time on these numerical integrations over a three-dimensional time-spatial super-surface as when we use standard numerical integration such as the Gaussian quadrature rule. So to avoid this numerical difficulty, we here adopt a kernel approximation [5] for the boundary integrals, by modifying the fundamental solution (9) in the part of the boundary integral domain where the condition I v f t ' / I r ' I << 1 is fully satisfied to

(9) &0 denotes the initial potential distribution prescribed a priori, t ' = ~"- t >/0. H stands for the Heaviside step function, imposing the causality. Following a similar process of discretization, boundary-solution approximation and boundary integration as in the steady state, we can also readily find the boundary element matrix equation from eq. (8) in a similar way as in the steadystate case [3,6].

The above, approximated kernel function can be analytically integrable with respect to time, because the expression is identical to the fundamental solution of the diffusion equation [14]. A numerical example is shown in Fig. 4 [5]. The boundary element solution (small circle) is compared with the exact solution (solid line). In spite of a relatively rough discretization, it is seen that reasonable solutions are obtained. However, as the P6clet number becomes larger, the solution becomes less accurate. Now we are going to find a more accurate and effective integration method. In boundary element analysis for unsteady problems, some alternative numerical approaches [15-17] without time-dependent fundamental solutions are proposed, such as the Laplace transform [15] method and the time difference method [17]. Although threedimensional analyses by means of these methods have never been presented, it is expected that they have a considerable potential in unsteady-state analysis of three-dimensional convective diffusion problems. In the future, the authors will investigate transient numerical solutions of the three-dimensional convective diffusion equation through numerical experiments employing these numerical approaches.

4.2. Boundary element solution In order to demonstrate the transient boundary element solution, as example, we shall treat the same domain as the steady-state problem shown in fig. 1. It is also assumed that the vector v = ( v x,O,O), v x= constant in time and space. When we solve the three-dimensional problem, we will readily that it is very difficult to analytically integrate the boundary integrals, because they include both

5. Conclusion In this paper, the author has presented three-dimensional boundary-element solutions for the convective diffusion equation with constant velocity. In particular, we have studied three-dimensional mixed boundary-element solutions. Further, it has been shown that both steady and unsteady solutions are stable and

T. Honma / The convective diffusion equation accurate up to comparatively high cell-P6clet numbers in comparison with the finite difference method and finite element method. Such a special convective diffusion equation with constant velocity could be applied in electrical engineering to traveling magnetic fields, such as eddy current problems in a linear induction motor [18]. In the near future, the authors will attempt to solve such field problems in electromagnetic solid mechanics [19] on the basis of the present work.

References [1] P.J. Roache, Computational Fluid Dynamics (Hermosa Pub., Albuquerque, 1976). [2] T. Ikeda, Maximum Principle in Finite Element Models for Convection-Diffusion Phenomena (Kinokuniya, Tokyo, 1983). [3] Y. Tanaka, T. Honma and I. Kaji, On mixed boundary element solutions of convection-diffusion problems in three dimensions, Appl. Math. Modelling 10 (1986) 170175. [4] Y. Tanaka, T. Honma and I. Kaji, On transient boundary-element solutions of unsteady convective diffusion equation, Trans. IEICE Jpn. E 70 (4) (1987) 242-244. [5] Y. Tanaka and T. Honma, An effective boundary-element analysis of three-dimensional unsteady convective diffusion equation using kernel approximation, Trans. IEICE Jpn. E 71 (4) (1988) 309-311. [6] M. Ikeuchi, M. Sakakihara and K. Onishi, Constant boundary element solution for steady-state convective diffusion equation in three dimensions, Trans. IECE Jpn. E 66 (1983) 373-376. [7] Y. Tanaka, T. Honma and I. Kaji, 0-1 order mixed boundary element solution in three dimensions, Trans. IECE Jpn. E 68 (1985) 409-410. [8] T. Honma, Y. Tanaka and I. Kaji, Regular boundary

319

element solutions to steady-state convective diffusion equations, Eng. Anal. 2 (2) (1985) 95-99. [9] Y. Tanaka and T. Honma, Boundary element solution of three-dimensional convective diffusion equation for high p6clet numbers, Proc. 1st Joint Japan/U.S. Symp on BEM, Oct. 3-6, 1988, to be presented. [10] M. Ikeuchi and K. Onishi, Boundary element solutions to steady convective diffusion equations, Appl. Math. Modelling 7 (1983) 115-118. [11] M. Ikeuchi, Transformed boundary element method for steady-state convective diffusion problem, Trans. IECE Jpn E 68 (9) (1985) 602-608. [12] Y. Tanaka, T. Honma and I. Kaji, Mixed boundary element solution for three-dimensional convection-diffusion problem with a velocity profile, Appl. Math. Modelling 11 (1987) 402-410. [13] M. Ikeuchi, Transient solution of convective diffusion problem by boundary element method, Trans. IECE Jpn. E 68 (7) (1985) 435-440. [14] Y. Tanaka and T. Honma, Three-dimensional boundaryelement method for magnetic diffusion model, Proc. 4th Joint MMM-Intermag Conf., July 12-15, BG-11 (1988), to be presented. [15] N. Okamoto, Analysis of unsteady convective diffusion problem by combined method with finite and boundary elements using laplace transformation, Proc. 1st Syrup. on Numerical Fluid Dynamics, Dec. 22-24, 1987, pp. 247-250. [16] M. Sakakihara, Characteristics incorporated boundary element for one-dimensional convective diffusion problem, in: Boundary Elements VIII, eds. M. Tanaka and C.A. Brebbia (Computational Mechanics Pub., 1986) pp. 793-801. [17] A. Taigbenu and J.A. Liggett, An integral solution for the diffusion-advection equation, Water Resour. Res. 22 (8) (1986) 1237-1246. [18] K. Ogawa and S. Nonaka, Boundary element analysis of single-sided linear induction motors, Trans. IEE Jpn. Sect. B 105 (12) (1985) 985-992 (in Japanese). [19] F.C. Moon, Magneto-Solid Mechanics (John Wiley & Sons, New York, 1984).