The auxiliary boundary element method for time dependent problems

The auxiliary boundary element method for time dependent problems

Journal of Computational North-Holland and Applied Mathematics 12&13 (1985) 239-245 239 The auxiliary boundary element method for time dependent ...

412KB Sizes 3 Downloads 28 Views

Journal of Computational North-Holland

and Applied

Mathematics

12&13 (1985) 239-245

239

The auxiliary boundary element method for time dependent problems G. De MEY Laboratory of Electronics, Ghent State University, 9000 Ghent, Belgium Received 4 June 1984 Revised 25 September 1984 Abstract: The boundary integral equation method or boundary element method is a well known technique to solve problems described by partial differential equations. By putting a source function on the boundary the partial differential equation is replaced by an equivalent integral equation defined on the boundary. Recently, some authors developped an alternative way by introducing an auxiliary boundary besides the physical boundary on which the source function is defined. Until now, this method has only been used for problems described by elliptic equations. In this contribution, it will be shown that the same idea can be used for parabolic equations. The technique has been applied to the diffusion equation. Experimental results have been obtained for a thermal diffusion problem including non linear boundary conditions. The experiments are also compared with the ‘classical’ boundary element method.

Keywords:

Boundary

element

method,

diffusion

equation.

1. Introduction During the last decennia, the boundary element method (BEM) has grown up to a powerful1 numerical technique, comparable to the finite element method or the finite difference approximation [6]. A lot of problems occuring in physics and engineering have been attacked using BEM [1,5,7,11]. The basic idea is using Green’s formula, capable to convert a problem described by the Laplace equation in an area S and appropriate boundary conditions on aS into an equivalent integral equation along as. Dividing the boundary into elements yields discrete versions of the integral equation. The same procedure can also be applied to time dependent problems described by parabolic equations. A fundamental feature of the BEM is that the discretisation takes place on as, i.e. the same curve (for a 2-D problem) where the boundary conditions are defined. One can imagine another numerical procedure with the discretisation on another curve, called the auxiliary boundary, inside or outside as. Using this idea some experiments on potential problems have been done by several authors independently from each other [4,10,13]. This explains the names “regular boundary element method”, “ auxiliary boundary element method”, “method of the fundamental solution”, . . . used for the same technique. Some theoretical investigations have also been done by Christiansen [3]. We will furtheron use the name “auxiliary boundary” due to a lack of any official nomenclature. 0377-0427/85/$3.30

0 1985, Elsevier Science Publishers

B.V. (North-Holland)

240

G. de Mey / Boundary element method

The reason to introduce an auxiliary boundary is quite obvious and mentioned by all the authors using it [4,10,13]. The BEM discretises on the physical boundary as, hence the largest numerical errors will occur on the boundary. It is also observed that the evaluation of the unknown function is more accurate for internal points of S. because the discretisation errors are smoothed out the larger the distance from the boundary. It is then quite obvious to move the discretisation (= boundary elements) away to an auxiliary boundary as, outside aS in order to get higher accuracies. This idea was confirmed experimentally. However a second phenomenon was observed: the larger the distance between &S and &S, the more the matrix turned out to be ill conditioned. In many case an optimum could be found [4]. As a rough approximation one can state that the optimum distance between aS and &S, is comparable to the length of the nearest boundary element. Nevertheless, it should be emphasized that the method is still in the experimental phase. The boundary element method has also been used to solve time dependent problems described by parabolic equations [2,14,15]. element method for a time In this paper, a first attempt to use the auxiliary boundary dependent diffusion problem will be described. More specifically a transient 1-D heat transfer problem with radiation boundary conditions will be investigated.

2. Fundamental equations and boundary conditions We consider an infinite plane with thickness 2 (Fig. 1). The initial temperature at t = 0 of the plate being To, the plate is cooled due to convection and radiation effects at the boundaries x= +l and x= - 1. The temperature T( x, t) satisfies the diffusion equation:

a2T/ax2

a*T

2=x

= air/at.

(1)

dT

! I i i

i j

I I

-a

,l

I I i i i

clxa

i Fig. 1. Plate thickness

2 immersed

in cooling

medium

with radiation

boundary

conditions.

G. de Mey / Boundary element method

On the boundary

241

x = 0 one has the condition:

aT/i3x = - BT - RT4,

(2)

- BT represent the cooling due to convection and - RT4 the radiation effect. The ambient temperature was taken T = 0. Due to the symmetry around x = 0 the problem can be simplified furtheron.

3. Integral equation In order to construct an integral equation function of the equation (1) [12] G(x,

r) =

replacing

(1) and (2), one needs the so called Green’s

~e-x2~4r 2fi

T(x,

The solution T(x,

t) is now written

as

t) = To +l(g(r’)[G(x

- 1, t - t’) + G(x + 1, t - t’)]dt’

(4)

0

where the symmetry on the boundary:

of the problem

zc=l = +g(t)

has been included.

+lg(t’)[

aG(x

+a;

From (4), the gradient

can be calculated

* - f’)).Y=,dr’.

The first term $g(t) is due to the singular behaviour of aG/ax on the boundary [8.9]. Note that g(t) is still an unknown 1-D function whereas the original problem is 2-D. Inserting (4) and (5) in the boundary conditions (2) yields aG(x

M> + @f’)( = -B

T,+

I -R

I

+

1, t -

t’)

ax

‘g(t’)[G(O,

)x_ldl’

t-t’)+G(2,

t-t’)]dt’

t-t’)+G(2,

t-t’)]dt’

J0 To+

‘g(t’)[G(O, J0

1 I

4.

(6)

(6) is a non linear Volterra integral equation in the unknown function g(t). Once g(t) has been determined, the temperature can be calculated by evaluating the integral (4). The integral equation (6) can be easily solved numerically by dividing the time axis into equal elements with length At (fj = iAr). In the interval (ti_*, ti), the unknown function is approximated by a constant gi. The integral equation (6) reduces to +gi+

h gjp(i-j)= j=l

-B

T,+

i j=l

gja(i-j)

1[ -R

To+

igja(i-j) j=l

1 (7) 4

242

G. de Mey / Boundary element method

A

0.06

1

Fig. 2. Temperature on the plate surface results obtained with the auxiliary BEM.

vs. time for At = 0.01; B =I

and R = (0.1, IO). The dashed

lines are the

where the coefficients (Yand /3 are given by cu(k) = ~;~+l’Ar[G(O, 7) + G(2, r)]dr

+ -j+

[ /@TiJGee-‘~~k”‘)Af

- &Ge-‘/kAf],

At time step ti, { gr . , . g,_i} are known from previous time steps and gi can be determined from the non linear equation (7), which has been solved by a Newton-Raphson iteration method. The temperature can be calculated by discretization of (4) in a similar way. Figure 2 shows the variation of the temperature T(1, t ) on the boundary vs. time for various values of At and R (full lines). The use of a piecewise constant approximation is not recommended as the only possible way. In some cases higher order approximations (piecewise linear e.g.) are preferred or even necessary to get sufficient accuracy [15].

G. de h4ey / Boundary element method

4. The auxiliary

243

boundary method

The discretisation procedure will no longer be done on the physical boundary x = _+1. The unknown function g(t) is now defined on the auxiliary boundary at x = 1 + a (Fig. 1). Until now the distance a can be chosen arbitrarely. The temperature is written as T(x, 2) = T, + @)[c(x

- 1 - a, t - t’) + G(x + 1 + a, r - t’)]dt’.

(10)

The temperature gradient is given by: aG(x - 1 - a, t - t’) + aG(x+l+u, ax ax

t-t’) (11)

Note that the above expression can be used for all internal points - 1 < x d 1, because the singularities of aG/ax occur only at x = f (1 + a). Problems related to singularities are eliminated by the auxiliary boundary. The construction of an integral equation and its numerical solution can be processed in the same way as has been done in the foregoing section. One has only to replace (4) and (5) by (10) and (11) respectively. Some results are shown on Fig. 2. The boundary temperature T(l, t) has been drawn vs. time for various values of At, R and the distance a (dashed lines). One observes that for every set of A T/l. tl

B= 1.

At=&01

l.-

0. 0.

Fig. 3. Deviation

I 0.1

I

0.2

*

t

A between auxiliary and classical BEM as a function of a/&.

244

G. de Mey / Boundary element method

parameters an optimal a-value can be found. If a is far from the optimal value large errors are observed especially at low times. This phenomenon can be understood physically. If CI is too high, the source function defined at a large distance a from the physical boundaries can hardly represent the exact behaviour of T(x, t). If the distance a is too small, the singular behaviour of G and aG/ax combined with discretisation errors become important. An obvious way to find the optimal distance a is to plot the deviation A between the auxiliary and classical BEM vs. u. The intersection of these plots with the u-axis yields the optimum value. However is a/& is ruled along the horizontal axis it is observed that all these plots intersect the horizontal axis around the same point a/& = 1 (Fig. 3). This result is not completely surprising because the quantity Ax2/At occurs by transforming the finite difference approximation in a dimensionless form. However, a theoretical proof has not been found yet. The method presented in this paper can be extended to two and three dimensional problems. The determination of the optimal distance a is then replaced by finding the optimal auxiliary boundary as has been done for static potential problems [4].

5. Conclusion In this paper the auxiliary boundary element method has been introduced for time dependent problems. It was found that the results coincide very well with those obtained by the classical technique. Experimentally a relation was found giving the optimum u-value as a function of the time step At.

References [l] C. Brebbia, The boundary element method for engineers (Pentech Press, London, 1978). [2] C. Brebbia, J. Telles and L.C. Wrobel, Boundary Element Techniques: Theory and Applications in Engineering (Springer, Berlin, 1983). [3] S. Christiansen, Condition number of matrices derived from two classes of integral equations, Math. Methods Appl. Sci. 3 (1981) 364-392. [4] G. Fairweather and R. Johnston, The method of fundamental solutions for problems in potential theory, in: C. Baker and G. Miller, Eds., Proceedings of the Durham Symposium on Treatment of Integral Equations by Numerical Method (Academic Press, London, 1982) 349-359. [S] D. Ingham, P. Heggs and M. Manzoor, Boundary integral equation solution of non linear plane potential problems, IMA J. Numer. Anal. 1 (1981) 415-426. [6] M. Jaswon and G. Symm, Integral equation methods in potential theory and electrostatics (Academic Press, New York, 1978). [7] G. De Mey, Integral equation for the potential distribution in a Hall generator, Electronics Letters 9 (1973) 264-266. [8] G. De Mey, An integral equation method for the numerical calculation of ion drift and diffusion in evaporated dielectrics, Computing 17 (1976) 169-176. [9] G. De Mey, Numerical solution of a drift-diffusion problem with special boundary conditions by integral equations, Comput. Phys. Commun. 13 (1977) 81-88. [lo] G. De Mey, Integral equations for potential problems with the source function not located on the boundary. Computers and Structures 8 (1978) 113-115. [ll] G. De Mey, Potential calculation in Hall plates, Adu. Electronics and Electron Phys. 61 (1983) 1-61. [12] P. Morse and H. Feshbach, Methods of Theoretical Physics; Part I (McGraw-Hill, New York, 1953) 191.

G. de Mey / Boundaqv element method [13] C. Patterson and M. Sheikh, A regular boundary element method Fluids 2 (1982) 239-251. [14] L.C. Wrobel and C. Brebbia, Time dependent potential problems, (Pentech Press, London, 1981) [15] L.C. Wrobel, A boundary element solution to Stefan’s problem, Elements, Hiroshima, 1983, 173-182.

for fluid flow,

245 Internal. J. Numer. Methodr

in: Progress in Boundary Elements Metho& Proceedings

5th Conference

on Boundary