SURFACE INTEGRALS APPROACH TO SOLUTION OF SOME FREE BOUNDARY PROBLEMS

ABSTRACT This paper is a continuation of the publication [1] where integral equation techniques were applied to the solution of a generalized Stefan problem. The regularization of the corresponding system of nonlinear integral Volterra equations offered here is quite different from that in [1], hence several new algorithms and numerical experiments. For consistency and easy reference we start this paper with sec.6.


INTRODUCTION
Boundary value problems for parabolic partial differential equations, in which one of the boundaries is not given, are generally known as moving boundary problems or Stefan problems.They appear in modeling fluid flow in porous media, diffusion with simultaneous chemical reaction, heat flow with phase transitions, etc.In treating these problems, one seeks both a solution to the differential equation and an unknown function which describes the position of the moving boundary.Two conditions are thus needed on the unknown moving boundary: one to determine the position of the boundary itself and the other to complete the solution of the differential equation.
This paper deals with the solution of the following one-dimensional (in terms of space variables) generalized Stefan problem: given the heat equation in a semi-infinite strip Ou o42u (0.1) c?---Ox 2 O, x>s(t), t>0; the initial and boundary conditions are (0.2) u(x, O) O, x > 0; (0.3) u(s (t), t) r(t), > 0; and the "extra condition" a measurement of the flux through the boundary is obtained from the experiment: (0.4) c) u (s (t) t) h(t), t>0, (s(0)=0), cX where r(t) and h(t) are given continuous functions.The simplest possible situation described by this problem is the melting of a semi-infinite solid.In this case u(x, t) represents the temperature distribution in the solid and s(t) represents the position of the moving solid-liquid interface.The problem consists of finding u(x, t) and s(t), given the initial temperature for all x >_ 0 and the temperature and flux on the moving boundary.The relation (0.4) replaces the heat balance requirement on the moving boundary which in the case of the classical one-phase problem with no source term in the boundary has the form: where X is a dimensionless latent heat of melting.
(i) The solution of one-phase problem with no source term on the boundary, given the experimental measurement (0.4), is trivial and we mention it here only for the record.
(ii) The two-phase problem, in which the temperature distributions in both the liquid (ul) and the solid (u 2 phases and the position of the moving phase have to be found, is given in the form: ' c9 X_, 1, 2 (x < S(t) for u 1, and x > s(t) for U2) > 0; U I(S (t), t) U 2 (S (t), t) O, > O; O3U 2 0U ]/2 -s '1 -s ddt S t>O, where i and Yi are non-dimensional parameters.In this case, if the (u phase is not accessible for experimental measurements, the problem may be put in the form (0.1) (0.4) by letting: where h(t) is given and the equations are solved for U2(X, t) and s(t).Once these are found, (0.5) may be solved for Ul, x (s (t), t).
(iii) There may be an unknown heat source or sink on the moving boundary.In this case (0.4) takes the form: (0.6) h(t) -(s(t), t) d s + q(t), t>O.
(iv) There is also a combination of (ii) and (iii), which gives a two-phase problem with a source in the boundary, still reducible to (0.1) (0.4).
(v) If the liquid phase moves with the velocity v(t), which may be a known or unknown function itself, a convection term must be introduced in the heat flow equation for the liquid phase: If the convection is caused only by different densities pl(Solid) and p2(liquid) across the phase-change boundary, then the conservation of mass for two incompressible phases gives us a condition on v (0.7) '1 dt ,o2v Later we shall see that whenever v(t) is given (either offered independently of the relation (0.7) or calculated from it no matter what the reason for the convection is), the problem (v) in any combination with the problems (ii)-(iv) is reducible to the problem described by (0.1) (0.4).In the case of (i) (v), v(t) may be considered unknown itself and the problem of finding it is again reducible to the framework of (0.1) (0.4).

SINGLE AND DOUBLE LAYER POTENTIALS: REDUCTION TO A SYSTEM OF INTEGRAL EQUATIONS
Based on the fundamental solution for a heat operator and Green's formula, there is an integral representation for the solution of the boundary value problem (0.1) (0.4), extended by 0 for < 0 and x < s(t) o To simplify the form of the domain of our problem we introduce the substitution y=x-s(t) in (1.1).We thus obtain anewdomain (x> 0, t> 0, wherey is again replaced by x and two types of surface integrals concentrated on the new boundary x O: (x+s(O-s()) ( 12 (') 4(t-'r) (1.3) V(x, t) = Jo ,] r e d'r, (X+S(t)-s('r)) (1.4) W(x,t)=lfoV('OIx+s(t)-s('C)-P('C)l 4(t--) v 2 "c) e henceforth referred to as single and double layer potentials respectively.Remarks.It should be noted that V(x, t) and W(x, t) with arbitrary, bounded, continuous density functions, exist and satisfy the heat equation inside the domain of our problem.The sum of these two potentials, however, does not satisfy the boundary conditions (0.3) and (0.4) without proper choice of the densities.Some properties of the surface integrals (1.3) and (1.4) that we use here, are discussed in Appendix 1.
We shall seek a solution to the problem (0.1) (0.4) in the form of a sum of these two potentials with as yet unknown density functions/.t(') and v('): (1.5) u(x, t) 2W(x, t) 2V(x, t), where V(x, t) and W(x, t) are the single and double layer potentials defined above; and for the unknown density in W we choose the function p('r) (see also (1.2), i.e. v p).
In [2] it was shown that V and W satisfy the zero initial condition, and for continuous density functions (1.5) is a classical solution to the problem (0.1) (0.4).
We need to find density functions such that expression (1.5) will satisfy the boundary conditions (0.3) and (0.4).
From (0.3), the continuity of V atx 0, and the jump relation for W (Appendix 1), we get or (1.6) This is a nonlinear Volterra integral equation of the second kind for the unknown function p(t).

REGULARIZATION
As it is proved in Appendix 2, W x is continuous at x 0 and the first integral in (1.9) is convergent.On the other hand the order of singularity of its kernel (1.10) is too high for this integral to be absolutely convergent.This fact makes the analysis of the system of integral equations (1.6) and (1.9) extremely difficult.In this paper and in the one that follows ([4]), we suggest two different regularization techniques for the integral operator in (1.9).
Since W (see (1.4)) satisfies Then, using (2.2) in conjunction with the boundary condition (0.4), we obtain the equation (1.9) now in a "regularized foma": As it is clearly seen, the regularization concerns only the first integral in (1.9).It also follows from the discussion above that for A --+ 0 + p(z) K(A, t, 7, p) dz-2A B(A, t, p) --+ p(z) K(O, t, 7, p) dz.0 Now, the system (1.6) and (2.4) remains to be solved.Remarks. 1) It is necessary to keep in mind that the value A, however small, must be chosen at the beginning of the regularization process and then it must remain permanently fixed.2) It is also worth mentioning that the system is nonlinear only with respect to one part of the solution, p(t), but this nonlinearity is extremely "heavy".(See also (1.2), (1.10), (2.3) and (3.3) below.)

EXISTENCE AND UNIQUENESS
We intend to establish that the matrix integral operator of the system (1.6) and (2.4) is a contraction in some closed ball U F from the space C(R +) (R) C(R+) of continuous vector-functions with the uniform vector norm (3.1) f = max {lip II, II/2 <F < o,,, where f= (p(t), la(t)) U F is a solution vector of the system.
We distinguish three integral operators in (3.2): (Due to (3.1), continuity of r, and an integrable singularity in (3.2), the solution of (3.2) also belongs to C(R+).Hence, A leaves the space C(R+) invariant).Studying these operators separately, we obtain the following estimates for any pair of vectors fl and f2 .t-y where kl, 2 k(t, v, Pl, 2 and II k k 2 II II Pl P2 II As a result, from the estimates above and the assumption (3.1), we finally get o< t<T Let us consider now the equation (2.4).The last integral in (2.4) can be expressed in the form 4 and the estimate similar to that of the operator A above can be obtained: The remaining integrals in (2.4) (see also (2.3)) represent a sum of space and time derivatives of the double layer potential W (1.4), that are calculated (due to our regularization technique) inside of the domain A > 0).This guarantees (e.g., see [2, 3]) their absolute convergence and continuous dependence of t, which in the linear case would immediately imply contractability of the corresponding integral operators.
To get the same results in our (heavily nonlinear) case several additional estimates for the singular integrals in (2.4) are required.For illustrative purposes we distinguish from (2.4) (see also (2.3)) three of the most singular cases: (+s(t)-s(z)) where A, as we mentioned above, is positive and fixed (analysis of I 2 is reducible to 13 case).
That is, for sufficiently small the right-side terms in (3.6) (3.7) can be made as small as necessary.The same is true for all remaining integrals in (2.4).
Using (3.5) and the techniques applied earlier, we can write an estimate similar to (3.4); and, finally, for the whole system (1.6) and (2.4) in the matrix form f= AId] + R, we get A 1 A If2 < Cf-f -f2 II, which implies that the operator A is a contraction in U F for < T and CT 1/2 < 1.

CONTINUATION
After the solution of the system (1.6) and (2.4) is found locally (for small), we can apply a continuation technique induced by [6].For illustration we consider only half of the system, the equation (1.6), and put it in the form: Let the pair (Pl(t),/.tl(t))satisfy (4.1) for all e [0, 0]; presumably there is another pair of continuous functions (pz(t), ]./2(t)) that satisfies the "shifted" equation P2(t) =f(t + ) + fo p2(7:) Kl[t-"c; s2(t $2( 27 , where s and s 2 are obtained from (1.2) for P and P2 respectively.Then it can be easily verified that the new pair of functions (P(t), M(t)) defined by the formula: P(t) =pl(t) for [O, 0] and=P2(t-oO for [o, 0 + ] (the same for M(t) and S(t) will satisfy the original equation (4.1) for all e [0, o + ].Of course, equation (4.2) has yet to be solved (again locally, for small ).The existence or non-existence of its unique solution will define our ability to expand the local solution of the system (1.6) and (2.4) into the new interval [a, o + b'].This procedure can be repeated until the solvability of (4.2) fails.We will address this question again in [4] in conjunction with numerical algorithms.

CONCLUSIONS
The problem with moving boundary (0.1) (0.4) is a generalization of several Stefan problems (i) (v) and their various combinations, and the solution of (0.1) (0.4) implies solutions to all those problems.
2) The function s(t) is presumed to be either continuously differentiable or H61der continuous with exponent oc > 1/2.By means of the single and double layer potentials (1.3) (1.4) (see also Appendices and 2), the solution of the problem in the form (1.5) is reduced to the solution of the system of two nonlinear Volterra integral equations (1.6) and (1.9) with respect to the unknown vector function f (p(t), /.t(t)).The function p is given by the formula (1.2), and f is presumed to belong to a closed ball U F in C(R+) (R) C(R+) with F 2G.
Applying a special regularization technique, we replace (1.9) by the equation (2.4) and subsequently prove that the matrix integral operator of the system (1.6) and (2.4) is a contraction on U F locally, for sufficiently small T, which implies the existence and uniqueness of the solution (in the class of functions and under the conditions specified above).The solution can be extended repeatedly into a larger interval by the procedure described in section 4, given that the "shifted" equation (4.2) can be solved.
In the paper [4] that follows we introduce yet another regularization trick and solve the corresponding system by several different numerical schemes.

APPENDIX
Keeping in mind that the following derivations are necessary in sections 3, we require that the function s(t) be either continuously differentiable or H/51der continuous with exponent o > 1/2 (which is equivalent to p('r) being continuous or H61der continuous with exponent o > 1/2).We shall temporarily keep the same densities as in (1.3) and (1.4) and discuss the properties of single and double layer potentials V and W in general terms.
1) First, we make a note that a single layer potential (1.3) has an integrable singularity; and in the case of a uniformly bounded density it converges absolutely and uniformly, which in particular implies its continuity at x 0.
2) Then we consider the behavior of the double layer potential W. For the reasons of technical simplicity we assume the density function v(v) is continuously differentiable, although being a H61der continuous of exponent o > 1/2 it is also sufficient.We will show that this potential has a jump discontinuity at x 0. Starting with (1.4) we define a new potential (?c, t): t [ (X + S(t) S('g) ) p('g') _ ] -(x+s(t)-s('O)2 (a.1) 7"//(X't)--2 fO 2(t-.)3/24_j e 4(,-) Now, performing the change of variables (a.2) y x + s(t)-s(v) we can express (a. 1) in the form: (a.3) (x, -x + s(t) e -y dy=----/-erf 2 for all x > 0; and for x convergent to 0 + (a.4) lim +(x, t) .vt) x-->O Using a similar technique based on (a.2) for x 0 and the observation that for the smooth function s(t), as was mentioned above, we obtain 7ff0, t)=_ u',,t____z.erf 2 s ( t ) ] The combination of (a.5) and (a.4) gives us a "jump formula" for the modified potential (c, t).
limo+ Y/(x, t)=----+ 7/(0, t) x---Now we use this result to prove that W(x, t) with variable density v('r) satisfies the same jump relation at x 0. We do this by showing that the difference W(x, t) fx, t) is continuous there.Using the assumption of differentiability of v(t) and the formula 1 (x+s(t)-s('O)2 .
2,/--ff 2, t-v j we observe that the difference W Z/'is equivalent to a single layer potential, continuous at x 0. This fact, combined with the existence of the two limits for g/'and W and the "jump formula," leads to 1 v(t), lira NTr, t)-N0, t)= lim +W(x, t)-W(0, t) =-x0 x--0 and finally to (a.6) lira +W(x, t) v(t) + W(O, t) x O -This property of the double-layer potential can be also obtained by the reference to the general theory of parabolic equations (e.g., see [2,3]), since the substitution of variables y x s(t) implies a new equation with one variable coefficient, namely (2.1).However, the general theory does not supply the formulas obtained above, which are essential for the regularization technique applied here and in the paper that follows.
Similarly, we can obtain the "jump formula" for V x" Both relations (a.6) and (a.7) play a key role in the reduction of the problem to the corresponding systems of integral equations.
APPENDIX 2 We need to prove that W x is continuous at x 0. (See (1.7) and (1.8) where this property was applied).As before, we consider '(x, t) and W(x, t) (defined in (a.1) and (1.4) respectively) with the restriction that v(t) and s(t) be continuously differentiable (or H61der continuous of exponent o > 1/2).
Using the mean value theorem for v(t) and s(t) (or their H/51der continuity), we can reduce the integral (a.8) to a sum of two proper integrals, three integrals of the "single layer potential type" (all of them are continuous at x 0) and one integral of the form: (x+s(t)-s(v)) z I C o v'(*) x2 e 4(t-'r) (t-.)3/2where 0 < v* < t.To study this integral we consider the following difference for > 0 and x close to 0: (x+s(t)-s('r)) 2 (a.9)I(x, t)-I(0, t) < C v'll 4,/--x ja ( t -x ' b ' ) 3/2 e 4(t-at) d'ro By adding and subtracting appropriate terms we can rewrite the integral on the fight-hand side of (a.9) as the sum of a double layer potential with constant density v(z) and a single layer potential V*.We denote this integral by U(x, t) and express it in the form: (a. 10) U(x, t)= x, t) + V*(x, t) =--erf (x+s(tl-s()) z 4(t-'r) 2p(v) p(v*) e dr.