Computatmnal Mechanlcs z8 (z996) 55-62 ~ Springer-Verlag ~996

Application of combined BEM-FEM algorithm in numerical modelling of diffusion problems E. Majchrzak

55 Abstract The paper contains the presentation of a certain algorithm which can be used in numerical modelling ofthermal diffusion problems - in particular in a case of nonhomogeneous domains. The method constitutes a composition of the finite element method (for the first domain) and the boundary element method (for the second domain). The way of these methods joining results from the heat flux continuity condition assumed on the contact surface. As will be shown, this approach essentially decreases the number of unknown parameters (the number of equations in final resolving system). The description of the algorithm is supplemented by the examples illustrating the possibilities of its application.

geometry of the mould and its thermal properties determine the course of solidification process. Presented in this paper approach allows to construct the numerical model which assures that the both domains/21 and B2 will be in exact way taken into account, at the same time the number of unknown parameters (temporary temperatures) is essentially less than in the case of 'homogenous' numerical method application. The method constitutes a composition of the FEM (for domain/21) and the 2nd scheme of the BEM (for domain/22). It should be pointed out that the problem concerning the domain/21 can be non-linear (thermophysical parameters of domain/21 can be temperature-dependent). 1 The way of FEM and BEM coupling results from the Introduction continuity condition (or conditions) given on the contact The thermal diffusion process proceeding in domain of surface. homogenous solid body is described by a partial differential In the paper the mathematical description of the problem equation (the Fourier equation) and adequate boundary-initial will be formulated, next the basic formulas concerning the FEM conditions. If we consider the non-homogenous objects (e.g. and the BEM algorithms for non-steady state diffusion problems casting solidifying in sand mould) then the mathematical model will be presented. In the final part of the paper the examples of thermal processes is more complicated. The system of partial of numerical computations will be shown. differential equations taust be taken into account, while the boundary condition on the contact surface results from the 2 continuity of normal heat fiux flowing through this surface. In Governing equations practice one can find the problems for which the domain/21 is A thermal diffusion process proceedings in heterogeneous more 'important' than domain/22 and the knowledge of domain/2 =/21 w/22 is described by a system of partial temperature fiel&, temperature gradients etc. in domain .621is differential equations of the form very essential, whereas the same information concerning the domain/22 is needless. On the other hand, however, the analysis X=.Q~: C OT~(X't)=div[2~grad Tl(X,t)] (1) of heat transfer processes in the system .r <;/22 requires the 1 Ot solution of the 'full' problem, because the peripheral domain /22 determines the thermal state of the 'main' domain/21. The typical example of such situation can be above mentioned C 2 a T z07 (X, t) = div[22 grad T2 (X, t)] (2) X = ~r system casting-mould. A designer of adequate technology is interested in details concerning the temporary temperature and the following boundary-initial conditions fiel& in casting volume, the kinetics of solidification process, the temporary shapes ofliquid metal sub-domain etc. From this point of view the mould domain is peripheral one, but the OT, (X, t) _ T1(X, t) - Y2(x, t) 8T2(X,t) On

R (X, t)

-

-

A2

OFl

(3)

Communicated by T. Cruse, 8 ]anuary 1996 E. Majchrzak Department of Mechanical Engineering, Silesian Technical University, Konarskiego 18a, 44-100 Gliwice,Poland

Correspondence to: E. Majchrzak This research work has been supported by KBN (Grant No 3 P404038 07).

Xel-'to:

[

(Pi T~,

OT~(X, t)] = O, ¦

[ aT~(X,t)] ~ j=o

xGG: % G t:o:

Z(x,o):%(x)

rdx, o): T~o(X)

(4) (s)

where C1, C2, 21, 96 2 are the thermophysical parameters (specific heats related to an unit of volume and thermal conductivities), ¦246 is a normal derivative at the point X e l ~, R(X,t) is a thermal resistance between sub-domains considered, F~0, F20 form the outer surface of the system. In the case R (X, t) = 0 the heat flux continuity condition resolves itself into the condition additionally determining a continuity of temperature fiel&

( ) ¦ Xe F12: 56

'I - - ~1~ ~ ~ - n

The mathematical model described in previous chapter in the further part of the paper will be assumed in the form

CI(T)

¦ T 1(X, t)

et

- V[21(T)VTI(X,t) ]

(7)

c OT2(X,t) 2 -~ -- 22V2T2(X,t)

(8)

,~2 ¦ =

c3n

(6)

pT1(X, t) = T2 (X, t)

X@F~2:

{q12(X, t) = - q21 (X,t) = q (X,t)

T2~(X,t)

(9)

T~~(X,t) - R ( X , t ) q ( X , t )

The domain considered is shown in Fig. 1. XeF~o: q~(X,t) = dL(X,t), As was mentioned, in a practice one can find the problems for which the domain f22 can be treated as a secondary XeF2o: q2(X,t) = ct2(X,t) (10) importance one (the details concerning the course of thermal processes in this domain are insignificant), but on the other hand the thermal interactions between ~~ and ~2 should be in t = 0 : TI(X,O)=Tm(X) T 2 ( X , 0 ) = 0 (11) exact way taken into account. Considering discussed problem it is possible (as a rule) to assume that the energy equation where q (X, t) = - fl 8 T (X, t)/3n. corresponding to ~2 sub-domain is a linear one (constant values It should be pointed out that the variables T1(X, t) and T~(X, t) of thermophysical parameters), while this simplification is not are transformed according to formula acceptable in the case of s domain. It should be pointed out Tk(X,t) --* Tk(X,t) -- T2(X,O), k = 1,2, at the same time it is that the numerical solution concerning O 2 will be found on assumed that initial temperature of ~2 domain is a constant the basis of the 2nd scheme of the BEM. This approach does one. In place of the Neumann conditions given on the outer not require the discretization of ~2 interior (if only the reference surface of the system (Eq. (10)) the other type of boundary level assuring zero initial condition will be assumed), and it is conditions can be also taken into account. a reason of essential reduction of unknown parameters (internal The numerical solution of discussed problem requires a time temperatures). mesh introducing 3 Numerical model In this paper the 2D model is considered. A heterogeneous domain ~ = g21w D~ is oriented in cartesian coordinate system X = {x,x2} and its discretization is shown in Fig. 2.

O=t~

l

(12)

and At = ts - t/ ~ is a step oftime. 3.1

The FEM algorithm for domain -Qt Let us consider the transition tf 1~ t{ Index s denotes a certain point Dom the interval At(e.g. t ' = t f y o r t ' = tf etc.), while w(X) = w(xl,Xa) is the so-called tapering function. The weighted residual method (Brebbia, Telles, Wrobel (1984)); (Mochnacki, Suchy (1995) used for Eq. (7) leads to the formula

{ V [ 2 1 ( T ) V T l ( X , t ) ] - C ~ ( T ) ~ ~cO~T-1;(X, , = vt)w~t

"X'; d,C21 = 0 (13)

Fig. l. A heterogeneous domain The last equation can be transformed to the form

[¦

~ t)TI(X,t) . Œ

~L--~Cx~ +

~ ¦

) ¦

~~ ¦

C1 cqTl(X't)l~-3 , , ' w(X)dS21= Yr, 96 1 ~

9w(X)dF~ Fig. 2. Discretization of domain 12

+-~~x~

where ~ = F n w/-]o"

1

ds

j,=,, ~

]t=t'

=

(14)

The s'2~ domain is divided into linear finite elements (triangles) - Fig. 2, and the Eq. (14) can be written in the form

8TI(X,t)] dOle ) f F8w(X)~ 8T~(X,t). ¦ E 3 ~1 -fr-A1 ~Xl ~X2 ~X2 Jt=t s (e) te) L ~Xl -

-

+ ~ f [ c , ST'(X't)~ (e) (e) L

~

w(X)d~~ y

Jt=t'

F ¦ = (b~)(~)L ~, ~

.Jt=pw(X)dF(lbe)

(15)

where I~(y is a square matrix (3 x 3) resulting from the geometry of finite element (coordinates of vertexes), while K(y 2(~~1/4[D I" ~(e)is a stiffness matrix (conductivity matrix (Mochnacki, Suchy (1995)); (Mochnacki, Suchy (1993)); (Tal-Ran Hsu (1994)) for finite element (e), 96 e) is a mean value of thermal conductivity, [Dr is the area of finite element, at the same time one can see that it was assumed t s = fr. The second integral:

F y

~~Ly

~

l~,~(~)d~~ ~~

57

where (be) are the sides of finite elements being in touch with other domain or environment. The temperature distribution T(y (X, t9 in sub-domain of finite element (e) results from the formula

XE(e):

T(y

N2

[(T,)~] N3]J(T,)~_|

[(r,);j

where N~, N 2, N3 are the linear shape functions (Tai-Ran Hsu, 1994), while the indexes 1,2,3 identify the vertexes of triangle considered. Analogously the time derivative for finite element (e) is expressed as 9

w

8tl, l X~(e):

8T~~)(X't')-[N]

N3] (8T]y

N2

(17)

\ 8t !~ !

8t

(21)

(16)

pY 8tAl

L\- X/

/3

where P(~~is a square matrix (3 x 3) resulting from the shape of finite element, P (y ]DIC(e)/129f' Ic)is a thermal capacity matrix and Cff) is a mean specific heat for element (e). It should be pointed out that the mean thermal conductivity and specific heat taust be calculated for each step of time, so the matrixes K(y p(e) are time-dependent. The last integrah

. F) 8TI(X,t) ]

and

' Ti -- "ri i\

---0.5L

'~t' ),

X.ff(e)" ~T~e)(x'ts)

N2

~[N 1

r N31\

(22)

/?2 /73] - 0.5 L m3J

(18)

_

at

= q~Y[fi,

At

I2

where ql b~) is a given boundary heat flux corresponding to segment (be) considered, while m, = 0 for 'internal' vertex and m, = 1, i = 1, 2, 3 for 'boundary' vertex, L is a distance between At //3 boundary vertexes. The last stage of algorithm consists in adequate coupling of The tapering function w is defined as follows (Mochnacki, Suchy matrixes K (e), p(y the matrixes (22) corresponding to (1995)); (Mochnacki, Suchy, (1993)) boundary conditions in matrixes Kf, pi and B for whole domain f2~ (Mochnacki, Suchy (1995)). In this way the Eq. (15) takes a form Xe(e): w(e)(x)=[/?~ /;2 173] N2 (19)

fl

LNd

(Kf+~tPf)'Tf=

AtlPf'Tf~ + 1 B.q~

(23)

The frst integral apearing in Eq. (15) is equal to

y F¦

(e)L

-

8T~(X,t) -

8X1

i -

-

Sw(X) +

~Xl

-

~(˜

=[/?1 /?2

-

~

1 ~:(y

]')3' 7 ~

-

8X2

21

aT~(X,t)] ~X2

at the same time the row matrixes of coefficients fl, appearing on the left hand and right-hand sides of resolving system can be reduced (Mochnacki, Suchy (1995)). The Eq. (23) will be also written in the form

d'Q~y

Jt=ts

[(T,)/] TVI

LIT:;! J

(20)

Wf.T(=

l ~ p f . T f - I + B.q f At '

(24)

The BEM algorithm for domain The basis of the BEM algorithm construction is the weighted residual criterion for Eq. (8) (Brebbia, Telles, Wrobel (1984)); (Majchrzak (1991))

¦

~
~t-

In the case of the second scheme of the BEM and zero initial condition, only the boundary F2 must be discretized. So, the boundary F 2 is divided into N linear boundary elements and integral Eq. (29) takes a form 1

N

t;

c(~')T2(~ ',tf) + ~ Z J J q2(X',t)T*(~',X',ff, t)dtdF2, ~2j l

1 (~,X, ff, t)d~2dt = O,

t o g22

1

t~

ff

=~2)~=,~j,~_ Tz(XJ, t)q*(g',XJ, tf t)dtdI'2j

5~2

612~--

(25)

C2 lf

58

,

T*({,X, ff, t) =4na2(ff _ t) exp

[

r2 1

4a2(t ~ -- t) '

2s=lJ~I

¦

B

(30)

Assuming that we consider the constant elements with respect to time, this means

te[f

T2(XJ,t) = T2(XS,t~), q2(XS,t)=q2(XJ, t ~)

~,t87

the time integration can be done in analytic way, namely

( ~,X, ff , t) ~

g'J---C22 y T*(g"XJ'ff't)dt=~Ei( r~ ~, dl 4~r2z \4a2At ]

(32)

l j q,(y X,,tf, t ) d t = l d,Jexp( zv = C2 te ~ 2~ r~

(33)

ax7

~X 2

]cosoq

ff , t) CO80{2 1

d

1

=8na2(t v - tyexp[

4a2(t F- t)

(27)

2

r~ ~ 4a2At ]'

t'

gf s = ~ t~lT* ( :',X',ff , t)dt 1 . - ts) ]}, -47c,~2{EZ~4az(t fr~t~x)]_Ei[_4a2(t~~ -

where (c.fi Fig. 3) (28)

d = (x 1 - y coscq + (x 2 -- ~2) cos%

1 tl

__ d,;

T~ ts_l)_]_ expl

r~

q2(X, t) T*( y

tF, t)dF2dt

where

to F2 t l:

=

TdX, t)q

(r

at the same time F 2 = ~20 w F w Coefficient c(~) = 1 for ~eg22, and c(~)e(O,1) for ~eG.

t~)_]}

(35)

1 tF

j" j

(34)

<-' = ~)lq*(~',x,,t~,t)dt

Using the properties of fundamental solution we obtain the following form of Eq. (25) (assuming zero initial condition for domain D-2) y171T 2 ( y

(31)

~n

I

=-<[

F2jH

- T2(X7, t)q* (~ ',X s, t f t)] d t d ~ ,

3T* (~,X, fr, t)

[¦

t~

(26)

r is a distance between point X = (x,,x 2) and point ~. = (~, ~2) in which the concentrated heat source is applied. Resulting from Eq. (26) normal heat flux is eqnal to

q*(~'X'tV't) = --22

lN

+ ~ ~ ~ ~ ~ [q2(XJ,t)T*(~',X;,ff, t)

where T* is the fundamental solution and it is a function of the form

(29)

Ei(.)

Ei(.)

=

is the exponential integral function of the form exp( - #) d~L

(36)

and it can be approximated by certain elementary function (Cody, Thacher (1968)). The final system of equations can be written in the form

N N c(~')T2(~J, ff) + 2 ~ g,jq2(Xqtf)dF~~j = 2 ~ z,jT2(XJ, tf)dF2j f-1

+ Z [z~-' r~(xJ, t') -gs s

~2 Fig. 3. Domain ~2

1

f = 1,2 ..... F, while index i idenfifies the successive boundary nodes.

(37)

~/P( v,,2, [TI'-'q ~ LT{~ 1J-- BI"q(

The integrals in above equation should be determined numerically (Bolkowski (1993)). The system of Eqs. (37) can be also written in the matrix form

~( s=l\

G . q f = Z.T~4- ~

Gi-'.q~

zf-'.T;-

=

)

2, Lq;lfl)

--G2"q2f+,=,~' [Zf-' Z(1-'] T 1

, f = 1,2 ..... F (38)

(43

3.3 The composition of the FEM and BEM algorithms For the needs of further considerations the following quantities will be distinguished (c.f. Fig. 1): - T (, q(denote the temperatures (boundary and internal) and heat fluxes on the outer surface of the ~x domain for t = fr, whereas T(-~-the same internal and boundar• temperatures for t = fr-l, - Tf, q[-denote the temperatures and heat fluxes on the outer surface ~0, - T~2, Tl:, q{2, q{1-are the temperatures and heat fluxes on the contact surface ~2" Non-steady temperature field in area D~ is searched on the basis of the FEM, and (for accepted denotations) the following resolving system is obtained (c.f. Eq. (24))

Iw( w(~]FT;] ~ [P( Pf] ['r;-'] 4- [B1 12] LT{d LT(?~J q(2

(39)

Non-steady temperature field in ~2 domain is found on the basis of the 2nd scheme of the BEM. So, for denotations introduced previously, the following system of algebraic equations should be taken into account (c.f. Eq. (38))

and T f, = T(2 -- qfR f

(44)

The solution ofmatrix Eq. (43) and formula (44) allow to determine the temperature fields in the D~ and ~Q2sub-domains for t = tf. 4 Examples of numerical simulation

As an example the numerical solution obtained for symmetrical fragment of rectangular iron casting (2D problem) with dimensions I0 x 50 [mm] solidifying in typical moulding sand will be presented. The following input data have been introduced: pouring temperature ofmolten metal Tl0 = 1395 ~ Iiquidus temperature TL = 1195 ~ (carbon equivalent: CE = 3.92%), border eutectic temperature T = 1145 ~ range of eutectic transformation AT = 5K, initial mould temperature T2o = 20 ~ and R = 0 (Majchrzak, Mendakiewicz (1995)). In numerical realization the temperature was defined as T ~ T -- T20.The substitute thermal capacity corresponding to finite element (e) has been defined as the mean value, at the same time T 1 < Teu - 0.5AT Teo - 0.5AT < T I < T~u + 0.5AT C,(T) = c3 TCu+ 0.SAT < T~ < TL Cl

[G:

Z21][T2f-[

Gj[q{]=[Z2

Lq~,

I C2

LT~d

C4

'-'

4-,~1([Z~'

Fl

z2f~']LTI, J

[G~'

G~I-' ,

[~713 L 21J/

(40)

The regard to continuity condition (9) transforms the above systems of equations to the following form

TI-_>rL

The values of cl, c2, c3 and c4 correspond to slopes of straight lines being the approximation of enthalpy-temperature diagram - Fig. 4 (resulting from enthalpy-carbon diagram

10.0

[W( W{2

[ Tlfl

1

-- B12] JT{2/ = ~ [P~

12 LT{~-lj +

Lqf•

B 1 "q(,

8.0

¦

IT{2] - % + z~~ w

4- 2

~=~\

[Z f-'

Z~?']

L'r{l

--IG[-'

G~~-s]

Lq;lJ/

(42)

W{2 --Z a

--B12

O ~ITi21

--G2, + Z21Rf - - Z J / q l [

L'l'~ l

3/ p:

~~I ~

E

Finally, the resolving system for the whole domain resulting from the coupling of (41) and (42) is of the form

[¦ If

~hase

6.0

- zd W / = - %ql

LT;~J

4

87

f PUll-Tl-'7/" /

(41)

[-z~l

(45)

"1- 4.0 aT~ 2.0

i

/

0j ~ 0

200

%i:: 400

600 800 T [~

1000 1200TL 1400

Fig. 4. Enthalpy-temperature diagram (CE = 3.92%)

59

N 6o iIXlNNN

Fig. 5. Discretization of domain f2 present in (Szargut (1977)), while the thermal conductivity 5~~(T) was assumed as constant value. The discretization of domain is shown in Fig. 5, while the results of numerical simulations are presented in Fig. 6, in particular the positions of distinguished isotherms for t = 1,3,5 and 10 [s] are marked. As the second example the numerical solution obtained for hot spot of T type will be presented. The following input data have been introduced: pouring temperature of molten meta] Tl0 = 1350 ~ liquidus temperature T87= 1300 ~ (carbon equivalent: CE = 3%), border eutectic temperature T u = 1145 ~ range of eutectic transformation AT = 5K, initial mould temperature T20= 20 ~ The discretization of domain is shown in Fig. 7. Additional set of points covering the interior of mould sub-domain allows to find the temporary temperature values in this domain, but as was mentioned this information

25

0

0

0

0

0

0

~

0

0

0

O

0

0

0

0

O

O

0 0

~ 0

0 0

0

0

0

0 0

O

0

0

0

O

~

O

0

0

0

0

0

0

~

Fig. 7. FEM-BEMgrid

is not necessary in order to determine the temperature field in casting volume. Figure 8 presents the positions of isotherms at selected times, namely t = 30 [s], 3 [min], 60 [min] and 70 [min].

t=30[s]

~,o¦

t/ltl l[[[d

Y&\%

9 tl llll I r

I~1

I

gl~ll

lillll~

!i ~ ~ ~~~

~ ~~-~~~~~~

IIIII1

lO i

III1!

;//~f 5

0

0

~ Alvt~

2345

Fig. 6. Position of isotherms (t = 1,3,5 and lO[s])

Fig. 8. Temporary position ofisotherms

-~

The correctness of presented algorithm was verified using Majchrzak, E. 1991: Application of the BEM in the Thermal Theory of Foundary, Publ. of the Sil. Techn. Univ., Mechanics 102, Gliwice experimental data quoted in (Mendakiewicz 1995). The Majchrzak, E.; Mendakiewicz, I. 1995: Numerical analysis of cast measurements of solidification process run concerned the iron solidification process, J. Mat. Process. Techn. 53, Elsevier: plates with different dimensions. The results of numerical 295-292 simulation (the first example) have been compared with the Mendakiewicz, I. 1995: Simulation of Cast Iron Solidification Process, experimental data (for the plate 10 x 50 [mm]) and the Doctor's Thesis, Gliwice compatibility from the practical point of view was quite Mochnacki, B.; Suchy, I. S. 1995: Numerical Methods in Computations satisfactory. So, it seems that proposed algorithm can be useful of Foundry Processes, Polish Foundrymen's Technical Association, for numerical modelling of heat transfer processes proceeding Cracow Mochnacki, B.; Suchy, J. S. 1993: Modelling and Simulation of Casting in non-homogeneous domains. Solidification, PWN, Warsaw Szargut, I. 1977: Numerical Methods in Thermal Calculations of References Industrial Furnaces, gl~sk, Katowice Bolkowski, S. 1993: Computer Methods of Electromagnetic Field Tai-Ran Hsu 1994: The Finite Element Method in Thermomechanics, Analysis, WNT, Warsaw Allen and Unwin, Boston, London, Sydney Brebbia, C. A.; Teiles, J. C. F.; Wrobel, L. C. 1984: Boundary Element Techniques, Berlin, New York, Tokyo, Springer-Verlag Cody, W. J.; Thacher, H. C. 1968: Rational Chebyshev approximations for exponential integral Ei(x), Mathematics J. Computations 22, 103: 641-649

61