The boundary element method in two-phase Stefan problems MARINELA TIBA Research Institute,for the Electrotechnical, 45-47, Tudor Vladimirescu Blvd, Bucharest, Romania INTRODUCTION In this paper we deal with the numerical treatment by the boundary element method of the problem: OH ----Ay =f in a T (1) 0t H(t, x) E [3(y(t, x)) in ~2 T (2) H(O,x) = Vo(X) in ~2 (3) Oy -- = u in 2; (4) On I2 is a bounded domain in R m , m t> 2 and I2 r is the cylinder I2x]0, T[ with lateral face ~=Px [0, T], P=0~2. We assume that voEL2(I2), u EL2(~), fE L2(I2T) and that/~ is a strongly monotone graph in R'R: (~y) -- ~(z))(y --z)/> t~ly --zl 2, ~ > 0 (5) bounded on bounded sets. When/~ is given by [ ~2(r--ro)+L, r>ro /~(r) =~ [0, LI, r = ro (6) / [al(r-ro), r<ro with cq, t~ 2 > 0, we obtain the two-phase Stefan problem) Many papers deal with the numerical solution of free boundary problems of Stefan type with Dirichlet condi- tions. 2-4 General boundary conditions are considered by White s, 6 by the finite differences method. Error estimates for the t'mite element method are obtained by Pawlow, 7 Jerome and Rose) The Neumann boundary condition was studied numerically and theoretically by the ffmite element method in Tiba and Tiba. 9 The boundary element method is based on the use of the fundamental solution of the corresponding equation. It is possible to define a fundamental solution for the non-linear problem (1)-(4), m but its form is compli- cated, with implicit character, and cannot be used. Therefore, we apply a special technique by using the fundamental solution of the Laplace operator. APPROXIMATION BY BOUNDARY ELEMENTS For numerical applications we consider [2 a polygonal domain in R 2. Then ~2 may be approximated by the union of the squares ~2k, k= 1,p, with disjoint interiors. Separately, we divide P in m segments F]. Let (Xi)i=l, m be Accepted April 1986. Discussion closes May 1987. the middle point of Fi and (X]+m)]=l, p be the gravity centre of I2/. We assume that the interval [0, T] is divided in N equal subintervals of length k > 0. We also assume that the data are given pointwise and we denote be if, u n the values off, respectively u, at moment n ,k and point x i. Let u~' be the fundamental solution of the Laplace equa- tion (Au* = --8), that is, in R2: 1 u*(x) = -- -- In Ix -- xi[, i = 1, m + p (7) 2n Taking Ui,i=l,m+ p as a test function and integrating twice by parts in equation (1), after an implicit discretiza- tion in time, we obtain the following integral relations: fu. -z n At u~ dI2 = [2 = __ C(Xi)yn+l(Xi) __ fyn+l OU* dr + d On F +fu*un+ldp+ffn+lu~dI2 (g) F ~2 Using the well known properties of the fundamental solution (7), it is possible to prove that c(xt)= 1 for interior points. When x i E P, c(xt) equals the angle made by P at x i (measured as units of 21r). Ifx t is a smooth point of F, then c(xt) = !2.n,12 We consider, that- the values of H n + 1, f n + x, y n + t are constant or each square I2k, k = l, p and given by the value in the centre of gravity. The values of yn and u n are supposed constant on each boundary element P/and equal with the value at the middle point of the element. For a given test function 'u*' and at the time moment 'n', equation (8)becomes: /-/k +m l-~{-t J u* da =- I2 k Ou* . . . . dP+ /=1 d On P/ f P[I"J~+m-t'¢n+l'~f'* + u'/÷' u* dr + e_,\ at -,k+,,,]j,,i d~2 1=1 k=l p/ ak (9) 0264-682X/87/010046-03 $2.00 46 Engineering Analysis, 1987, Vol. 4, No. 1 © 1987 Computational Mechanics Publications