- Email: [email protected]
A mathematical model for the freezing of a water-saturated porous medium
91
In computation (c) models Q see [l] and P were used. In Figs. 3,a and 3,b, respectively we show for these models the fields of r. v, p, V (V is the velocity modulus) in the meridional plane of flow (x, r) with x axis of symmetry at the instant t = 1.33. We measured Tin eV, v in c&/g, p in lOi g/cm.se?, V in km/set, and x and r in cm. All the calculations showed that the use of the law of conservation of energy (1) with the “isothermal jump” and “Knudsen layer” models leads to slightly different results. At the same time, the use of approximate energy equation (2) in both models of vapourization kinetics gives a difference of roughly twice in the vapourizing mass of aluminium as compared with calculations based on Eq. (1). The qualitative difference between tbc models shows itself, e.g. in models P and Q, as distinct from R and S, admitting vapourimtion due to interaction of the vapour and the surface in accordance with (3) and (4), even when there is over-all heat removal (F + fGTdi@ < 0).
Trmshd
by D. E. B.
REFERENCES 1. 2.
3.
ZUBOV, V. I., ef uf., Computation of laser action on a plane barrier and its vapour, Zh vychid. Mat. mpf. Fiz.. 23,6,1X&1522, 1983. MAZHUKIN, V. I., and PESIXYAKOVA, 0. A., Numerical modelling of processes of metal surface vapourization by laser radiation, Preprint IPMatem. Akad. Nauk, No. 48, Moscow, 19W. An algorithm for the numerical solution of the problem of surface vapourization of material by laser radiition, Zh. vychid Mat. mar. Fir, 25, 11, 1697-1709,198X BERGEL’SON, V. I. CIal., Formation of plasma in a layer of vapour arisii from the action of lasar radiation on solid, Kvanrovayu Elcktronika, l&4,20-27,1913.
4. 5.
KNIGHT, C. J., Theoretical modellmg of rapid surface vaporization with back pressure, AIM Jow~l, 17,5,519-523,1979. ZUROV, V. N. et al.. Calculation of the motion of the vapour of a solid under the action.of laser tad&ion, in: Rdiaring Gas Dynamics (Dinamika izluchayushchego gpzp), No. 3 Ws AN SSSR, Moscow, 1980.
U.S.S.R. ComPut. Maths. Matk Phys.. VoL 26, No. 6, pp. 91-95.1986
0041~5553/86 510.00+0.00 0 1988 Pergamon Journals Ltd.
Printed in Great Britain
A MATHEMATICAL
MODEL FOR THE FREEZING OF A WATER-SATURATED POROUS MEDIUM* A. M. MAKSIMOV and G. G. TSVPIUN
The problem of the freezing of a water-saturated porous medium is considered. It is shown that using the assumptions of the existence of a phase transition front and of the dependence of the temperature of crystallization of the water on the pressure results in supercooling in some domain of parameters. To construct a solution that does not allow su~r~oling, a two-phase zone model is introduced. The solution of many engineering geocryology problems is based on modelling the freezing process of the earth. To describe heat processes, as a rule, a well-known formulation of Stefan’s problem is used, that assumes the existence of a phase transition from [ 11. The propagation of a plane front of the freezing of a nondeformed water-saturated porous medium leads to water being forced out in the direction of motion of the front following an increase in volume according to the phase transition. Such filtration flow is described by Darcy’s law and consequently, assuming water is weakly compressible, by the piezoconductivity equation [2]. In a porous medium of low permeability the motion of the crystallization front will be accompanied by a significant increase in the pressure of the liquid. In this case the dependence of the local crystallization temperature on the pressure follows naturally. In the one-dimensional formulation the basic equations are the equations of heat conduction for a solid phase zone (the ground is completely frozen, o
and the equations of heat conduction and piez~o~ducti~ty ~T/&=e,iFTlW,
(1)
for the liquid phase zones (thawed ground, X(t)=GS-) w,i&-%~P~d?9.
(2)
Here x is the coordinate, t the time, T the tempeature, X(t) the coordinate of the phase transition front, P the pressure, a, and ar are the coefftcients of thermal diffusivity for the frozen and thawed ground respectively, x=k,K/tm is the coefficient of piezo-conductivity, k, the permeability of the thawed ground, K the bulk modulus of water, WIthe porosity, and p the viscosity of water. Note that as the values of the characteristic quantities show, convective heat transfer can be disregarded. The boundary conditions are integral analogues of the equations of conservation of heat and mass of the water, described for a moving interface of the zones X(r): T_-T+-T.-j(P.-Pe),
lZh. vych&l. Mat. mat. Fir.. 26, II, 1743-1747, t986.
(3.4
92
t3.W and for x=0 T-TO-q(O).
(4)
Herefis the linear dependence of the temperature of the phase transition on the pressure: f=T,-t&P.-P,,), where n is the coefficient of drop in temperature of the phase transition, T, is the crystallization temperature of water not under any excess pressure, A, and A, are the thermal conductivity of the frozen and thawed ground respectively, pi and b are the density of ice and water respectively, and the indices +, - and * denote values on the right, on the left and on the front respectively. The initial conditions are: TN,
0)-d',,
P(t, O)-PO,
TOW(O),
(5)
X(0)-0.
If P0 is a constant, then without loss of generality we can put PO=0 and use the quantity P to mean excess pressure that causes a drop in temperature for water crystallization. If in the initial and boundary conditions (4). (5) the quantities To and 10 are constant, then the problem has the self-similar solution r-T(k),
P-P(E),
f -rl-‘l*.
X-+‘~,
(6)
In the solid phase domain @C&9) this solution has the form T-T”+
(T.-Ta)erf(~/2a,“)erf-‘(S/2a.“l),
(7)
and in the liquid phase domain (@<&GQ) T-T,+
(T.-T,)erfc(E/Zal’~)erfc-*(B/Za,”),
P-P.
(84
erfc(E/Zx”)erfc-1(B/2x’“).
VW
After substituting (t?)-(8) into (3) we have a system of equations that defines the quantities T., P. and 8:
By eliminating T. and P. from this system we obtain a transcendental numerically.
FIG. 1
equation for & which can easily be solved
FIG. 2
Figure 1 represents the dimensionless temperature p=T-T,)/6 (the solid line), pressure &P/B (the dashed line) and the local temperature of the phase transition %-qPA3 (the dot-dash line) as a function of the dimensionless selfsimilar variable ~=~/o,r, worked out using the following values of the parameters for the ground [3]: T0=274%, TO=263%, o,=0.78X10dmz/s, u,~l.14X10dm2/s, x=0.86X10dm2/s, mz0.2. Here &q&h, is the characteristic temperature, where u, and A, are the thermal diffusivity and thermal conductivity of water respectively. In all calculations of the parameters that define the properties of the water and ice, the following values are used 141:a,=2.3XIOdm2/s, A,=O.S8W/(m.*K), C,,=4.19kJ/kg.°K, &2XfO*Pa, ~=O.OOlPa.s, n=0.765XlO~‘%/Pa, z=273”K, q=333.7kJ/kg, p~=1000kglm3, p,=910kg/mf, C,=2.09kJ/kg.°K, A,=2.23W/m.OK. Here C, and C, are the heat capacity of water and ice respectively, A, is the thermal conducti~ty of ice, and q is the heat of the phase transition. It is characteristic that the curve of the self-similar solution Tin the liquid phase zone everywhere lies above the curve of the local phase transition temperature Tp However a domain of values of the parameters exists for which the
93
dependence of the crystallization temperature on the pressure leads to supercooling of the fluid. An example of the calculation of such a solution is represented in Fig. 2 and corresponds to the following values for the ground parameters: T0=274W, P=263%, a,=0.78X10dm2/s, a=1.14X10d/s, x=0.86X10~‘m2/s, m=0.2 (the meaning of the curves is as in Fig. 1). As is well-known, in an analogous supercooling a frontal solution of the problem of crystallization of a binary melt is also necessary [S-7]. This problem is solvable if it is assumed that between the solid phase zone and the liquid phase zone there is a two-phase state zone in which water and ice exist in a state of thermodynamic equilibrium. In this case the local temperatures of the phases correspond and are equal to the crystallization temperature for local pressure values. The mathematical formulation of the problem has the form (1) in the solid phase domain (ocr
“WE
kc,
ap
LI
al
---(
T-f(P)-
Tc-qp,
P”“P2(lC$).
(9c)
Here pm,C, and A, are the density, heat capacity and thermai conductivity in the two-phase zone respectively, Y is the humidity (the volume portion of the Iiquid phase per unit volume of the two-part space), v. is the rate of filtration of the water, I& is the permeability in the two-phase zone, and pWois the density of the water when there is no excess pressure. The permeability of a porous medium is related to the humidity v by the following relation that is an approximation from experimental data [3]:
Moreover, it is assumed that the thermal capacity and the thermal conductivity in a two-phase zone obey the mixing rule: plllC~-(i-m)p,SC.~+mvp~~+m(l-v)p~i, L-(l-m)h.,+mvb,+mji-v)t
ppt C, and & are the density, heat capacity and therma conductivity of the surface of the ground respectively. The boundary conditions are integral analogues of the equations of conservation of energy and mass of the water, recorded for the boundaries of the zones taking into account that at the boundary between the liquid phase zone and the two-phase zone L’s(t) the humidity Y returns to I and at the boundary between the solid phase and the two-phase zone Z,(t) it takes the required value Y.:
Here
r-Xl(t):
v-v.,
T_-T+-T.-To-tjP.,
The initial conditions take the form T@, 0) -To,
P(z, O)=O,
T&f(O)-Th
x, (0)-X*(O) -0.
(11)
If in the boundary and initial conditions (10) and (11) the quantities To and 7@are constant then the problem has a self-similar solution of the ordinary type: T-T(E),
P-WE), x, (1)+‘I,
v-v(E), XI(i) -IW.
+rr-‘h,
(W (12b)
This reduces to solving a boundary value problem for the system of ordinary differential equations that have been set out clearly in (7) for the solid phase domain (KKy) and in (8) for the liquid phase domain f&C%=) for Tr=T,-#.
94
The constants y, /3, P’ and T. have to be found from the conjugation condition of the solution through the two-phase state zone in which after substituting (12) into (9) we have h,T”+T’v’(h,-hr)nr + $
IT’p,C,E+mTv’(pC,-p,C~)E+p,m~v’El-
0,
(13a)
At the boundary [=y the conditions I-T.-T,-qP.,
P’=
\-
A,T_‘-A,(v.)T+’
- +
mt~p,r.y,
are satisfied and at the boundary [=& the conditions T-T’-T,-qP*,
!!=i,
T-‘-T+‘,
P-‘-P+‘.
(15)
hold. The non-linear boundary-value problem (13x15) is solved numerically. In Fig. 3 the results of the calculation of the distribution T(the solid line), P (the dashed line) and Y(the dot-dash line) are given for the following values of the soii parameters [4]: Z&274%, ?%?63%, a,=0.7gX10dmz/s., u,~l.14X10dm2/s, ~=0.86X10“mz/s, nH.2, k,=8.SX10‘“mz, pS=2000kg/m~, C,=0.92kJ/kg°K, A,l=2W/m°K.
\f ._.-----. 6 Flu2v
\
.A-\
g
1 ‘1, ‘.
1s
6 06 4
.4
FIG. 3
The results show that the frontal model usuahy used for engineering caIcuiations can lead to an inadequate description of the freezing process of a material of low permeability. The realization of some or other freexing state depends qually on the value of the Permeability and on the values of the initial and boundary temperatures. These Parameters can be represented in nondimensional form as follows:
FIG. 4 In the plane of the Parameters (VO,& (Fig. 4) the position of the boundaries of the domain of realization of the freezing regimes for diierent values of $ is shown. The domains lying below the curves correspond to solutions with a phase transition front and those above the curves correspond to values of the parameters for which a two-phase zone is created. A decrease in $ leads, of course, to a widening of the domain of existence of the two-phase zone. Trandared by S. R.
I
95
REFERENCES 1.
TIKHONOV A. N. and SAMARSKII
A. A., 77~ Equations o~~azke~ticu~Phy~ic~
(Uravnenie ~temati~b~koi
fiziki), Nauka, Moscow,
1966.
2. 3.
4. 5. 6. 7.
COLLINS R, Flow ofF&d Through Porous Muteria~s, Mir, Moscow, 1964. Thermophysical Research on the Cry&e Zone of Siberia (Teplotizicheskic issledovaniya kriolitozony Sibiri, Nauka, Novosibirsk,1983. KUKHLING KH. Han&ok ofphysics, Mir, Moscow, 1985. IVANTSOV G. P., “Diffusion” supercooling in the crystallization of a binary alloy, Dokl. Akud Nauk SSSR, 81,2, 179-181.1951. BORISOV V. T.. The crystallization of a binary alloy with conservation of stability, Dokl. Akarl Nauk SSSR. l&3,583-586, 1961. ABDONIN N. A., A ~athe~ticuI Description of the Processes of Crysta#ization (Matematiche~oe opisanic protsessov kristaiiiaatsii), Zinatne, Riga. 1980.
(NW-5553/86 glO.~.~ 0 1988 Pergamon JournalsLtd.
U.S.S.R. Compztt. Maths. Moth. Phys., Vol. 26, No. 6, pp. 95-98, 1986. Printed in Great Britain
EXACT ESTIMATES OF THE ERROR GRADIENT FOR PROJECTION-DIFFERENCE FOR PARABOLIC EQUATIONS IN AN ARBITRARY DOMAIN*
SCHEMES
I. D. TURETAEV Estimates are given for the rate of convergence of some two-layer economic projection~ifferen~ schemes @.d.s.) and Crank-Nicholson p.d.s. for solving parabolic equations in an arbitrary domain. 1. An optimal estimate of order O(r + AZ)is obtained in [1] in the norm of h(Q) for the error of a family of economic p.d.s. for solving parabolic equations in an arbitrary domain fl c R*. In the present paper we prove for these p.d.s. orderwise exact estimates of the rate of convergence O(s + A) in the norm II,II(*J-IID~.R~+IID.u*. _. In addition, we study the Crank-Nicholson p,d.s. in an arbitrary domain fl c I?“, n Z 1. For the error of this scheme we prove exact estimates order order O(# + ]k]*) and O(S + ]hl) in the norms of f&Q) and ll.@r’respectively. For the case of constant coefftcients, when ff is a parallelelpi~, the optimai estimate of the error of the Crank-Nicholson scheme in the norm of L,(Q) was obtained in [2]. The rate of convergence of the p.d.s. whenfsatisfies some smoothness condition was studied in [3,4] (see also the references there). We solve the first initial value problem Dlu+Lu-f
in
Q-$2X(0, T),
ulho-uo(+),
u/r-o,
I%%aX(O. T).
(1)
where ~~-D~(Q~~(z)D~~) -+I u,r=Q (summation from f to n is understood with respect to repeated subscripts i, j), and fl is an arbitrary bounded domain of R: n L 1, with dfIcC. We assume that ~l~i*catrE,E,~r~lll* for any f = (51, * . . . 5.) = R, &%I (y&r, 1.1is the eikonal of the norm in R”, and o,~-o~(, an=W-‘(Q), o=&.(R), 6 I -4, 2,. . , n. Under these conditions we have for the solution of problem (1), (2) the estimate llD,Pull,+liD*Dtull~+llLII),nli~, ,
Here, Drp is the gradient of 9, D29 is a matrix with
D,Dw, ndil~-HD~f~l*+iif~n, +ILc& elements (here, are the norms in L*(Q) and f,Dcf=&(Q), fo, .a,&=U+(Q); Il.ilp. 11.11, @*(n) respectively), and K is a constant, independent of the functions in the inequality (and of r, A, if the data of the problem depend on them). We assume below in this section that n = 2. We construct a mesh domain (see [ 1, 51). Let R? be a rectangular mesh in R2 with steps A,, !tz with respect to the variables x,, x, respectively. Let &, be the domain which is the union of partial rectangles of Rh2which lie in fi and are spaced 0 (max (h2, !I,)) from the boundary XL Further, the domain a\& is divided into triangular elements in the usual way. Then, the required mesh domain f& is defined as the union of all triangles and f&. We assume that the division of fl into finite elements is regular. With this way of constructing fh, the boundary dfl is approximated by the boundary dfl, of the domain f& with order O(K), where A is the maximum of the sides of the elements. Let s” be the space of functions, continuous in fi, linear in triangles, and bilinear in the rectangular elements, and vanishing in ii/n,; S, is the space of functions, continuous in [O, r], linear in each interval fr(k-tf. rk], k=i, 2.. . . ,.4&TM-T. We put S = s” X S. We consider the following family of schemes for the approximate soiution of problem (1): (yr”‘, cp)+‘WDy?“, +a‘acym,
‘Zh. vychisl. Mat. mat. Fir., 26,
11,1748-1751,
1986.
Dcp)+O+(&Dxytm,
DiDtcp)p~
cp)==(~fl”. CPL m-o, 2, . . . ,Y-i,
In computation (c) models Q see [l] and P were used. In Figs. 3,a and 3,b, respectively we show for these models the fields of r. v, p, V (V is the velocity modulus) in the meridional plane of flow (x, r) with x axis of symmetry at the instant t = 1.33. We measured Tin eV, v in c&/g, p in lOi g/cm.se?, V in km/set, and x and r in cm. All the calculations showed that the use of the law of conservation of energy (1) with the “isothermal jump” and “Knudsen layer” models leads to slightly different results. At the same time, the use of approximate energy equation (2) in both models of vapourization kinetics gives a difference of roughly twice in the vapourizing mass of aluminium as compared with calculations based on Eq. (1). The qualitative difference between tbc models shows itself, e.g. in models P and Q, as distinct from R and S, admitting vapourimtion due to interaction of the vapour and the surface in accordance with (3) and (4), even when there is over-all heat removal (F + fGTdi@ < 0).
Trmshd
by D. E. B.
REFERENCES 1. 2.
3.
ZUBOV, V. I., ef uf., Computation of laser action on a plane barrier and its vapour, Zh vychid. Mat. mpf. Fiz.. 23,6,1X&1522, 1983. MAZHUKIN, V. I., and PESIXYAKOVA, 0. A., Numerical modelling of processes of metal surface vapourization by laser radiation, Preprint IPMatem. Akad. Nauk, No. 48, Moscow, 19W. An algorithm for the numerical solution of the problem of surface vapourization of material by laser radiition, Zh. vychid Mat. mar. Fir, 25, 11, 1697-1709,198X BERGEL’SON, V. I. CIal., Formation of plasma in a layer of vapour arisii from the action of lasar radiation on solid, Kvanrovayu Elcktronika, l&4,20-27,1913.
4. 5.
KNIGHT, C. J., Theoretical modellmg of rapid surface vaporization with back pressure, AIM Jow~l, 17,5,519-523,1979. ZUROV, V. N. et al.. Calculation of the motion of the vapour of a solid under the action.of laser tad&ion, in: Rdiaring Gas Dynamics (Dinamika izluchayushchego gpzp), No. 3 Ws AN SSSR, Moscow, 1980.
U.S.S.R. ComPut. Maths. Matk Phys.. VoL 26, No. 6, pp. 91-95.1986
0041~5553/86 510.00+0.00 0 1988 Pergamon Journals Ltd.
Printed in Great Britain
A MATHEMATICAL
MODEL FOR THE FREEZING OF A WATER-SATURATED POROUS MEDIUM* A. M. MAKSIMOV and G. G. TSVPIUN
The problem of the freezing of a water-saturated porous medium is considered. It is shown that using the assumptions of the existence of a phase transition front and of the dependence of the temperature of crystallization of the water on the pressure results in supercooling in some domain of parameters. To construct a solution that does not allow su~r~oling, a two-phase zone model is introduced. The solution of many engineering geocryology problems is based on modelling the freezing process of the earth. To describe heat processes, as a rule, a well-known formulation of Stefan’s problem is used, that assumes the existence of a phase transition from [ 11. The propagation of a plane front of the freezing of a nondeformed water-saturated porous medium leads to water being forced out in the direction of motion of the front following an increase in volume according to the phase transition. Such filtration flow is described by Darcy’s law and consequently, assuming water is weakly compressible, by the piezoconductivity equation [2]. In a porous medium of low permeability the motion of the crystallization front will be accompanied by a significant increase in the pressure of the liquid. In this case the dependence of the local crystallization temperature on the pressure follows naturally. In the one-dimensional formulation the basic equations are the equations of heat conduction for a solid phase zone (the ground is completely frozen, o
and the equations of heat conduction and piez~o~ducti~ty ~T/&=e,iFTlW,
(1)
for the liquid phase zones (thawed ground, X(t)=GS-) w,i&-%~P~d?9.
(2)
Here x is the coordinate, t the time, T the tempeature, X(t) the coordinate of the phase transition front, P the pressure, a, and ar are the coefftcients of thermal diffusivity for the frozen and thawed ground respectively, x=k,K/tm is the coefficient of piezo-conductivity, k, the permeability of the thawed ground, K the bulk modulus of water, WIthe porosity, and p the viscosity of water. Note that as the values of the characteristic quantities show, convective heat transfer can be disregarded. The boundary conditions are integral analogues of the equations of conservation of heat and mass of the water, described for a moving interface of the zones X(r): T_-T+-T.-j(P.-Pe),
lZh. vych&l. Mat. mat. Fir.. 26, II, 1743-1747, t986.
(3.4
92
t3.W and for x=0 T-TO-q(O).
(4)
Herefis the linear dependence of the temperature of the phase transition on the pressure: f=T,-t&P.-P,,), where n is the coefficient of drop in temperature of the phase transition, T, is the crystallization temperature of water not under any excess pressure, A, and A, are the thermal conductivity of the frozen and thawed ground respectively, pi and b are the density of ice and water respectively, and the indices +, - and * denote values on the right, on the left and on the front respectively. The initial conditions are: TN,
0)-d',,
P(t, O)-PO,
TOW(O),
(5)
X(0)-0.
If P0 is a constant, then without loss of generality we can put PO=0 and use the quantity P to mean excess pressure that causes a drop in temperature for water crystallization. If in the initial and boundary conditions (4). (5) the quantities To and 10 are constant, then the problem has the self-similar solution r-T(k),
P-P(E),
f -rl-‘l*.
X-+‘~,
(6)
In the solid phase domain @C&9) this solution has the form T-T”+
(T.-Ta)erf(~/2a,“)erf-‘(S/2a.“l),
(7)
and in the liquid phase domain (@<&GQ) T-T,+
(T.-T,)erfc(E/Zal’~)erfc-*(B/Za,”),
P-P.
(84
erfc(E/Zx”)erfc-1(B/2x’“).
VW
After substituting (t?)-(8) into (3) we have a system of equations that defines the quantities T., P. and 8:
By eliminating T. and P. from this system we obtain a transcendental numerically.
FIG. 1
equation for & which can easily be solved
FIG. 2
Figure 1 represents the dimensionless temperature p=T-T,)/6 (the solid line), pressure &P/B (the dashed line) and the local temperature of the phase transition %-qPA3 (the dot-dash line) as a function of the dimensionless selfsimilar variable ~=~/o,r, worked out using the following values of the parameters for the ground [3]: T0=274%, TO=263%, o,=0.78X10dmz/s, u,~l.14X10dm2/s, x=0.86X10dm2/s, mz0.2. Here &q&h, is the characteristic temperature, where u, and A, are the thermal diffusivity and thermal conductivity of water respectively. In all calculations of the parameters that define the properties of the water and ice, the following values are used 141:a,=2.3XIOdm2/s, A,=O.S8W/(m.*K), C,,=4.19kJ/kg.°K, &2XfO*Pa, ~=O.OOlPa.s, n=0.765XlO~‘%/Pa, z=273”K, q=333.7kJ/kg, p~=1000kglm3, p,=910kg/mf, C,=2.09kJ/kg.°K, A,=2.23W/m.OK. Here C, and C, are the heat capacity of water and ice respectively, A, is the thermal conducti~ty of ice, and q is the heat of the phase transition. It is characteristic that the curve of the self-similar solution Tin the liquid phase zone everywhere lies above the curve of the local phase transition temperature Tp However a domain of values of the parameters exists for which the
93
dependence of the crystallization temperature on the pressure leads to supercooling of the fluid. An example of the calculation of such a solution is represented in Fig. 2 and corresponds to the following values for the ground parameters: T0=274W, P=263%, a,=0.78X10dm2/s, a=1.14X10d/s, x=0.86X10~‘m2/s, m=0.2 (the meaning of the curves is as in Fig. 1). As is well-known, in an analogous supercooling a frontal solution of the problem of crystallization of a binary melt is also necessary [S-7]. This problem is solvable if it is assumed that between the solid phase zone and the liquid phase zone there is a two-phase state zone in which water and ice exist in a state of thermodynamic equilibrium. In this case the local temperatures of the phases correspond and are equal to the crystallization temperature for local pressure values. The mathematical formulation of the problem has the form (1) in the solid phase domain (ocr
“WE
kc,
ap
LI
al
---(
T-f(P)-
Tc-qp,
P”“P2(lC$).
(9c)
Here pm,C, and A, are the density, heat capacity and thermai conductivity in the two-phase zone respectively, Y is the humidity (the volume portion of the Iiquid phase per unit volume of the two-part space), v. is the rate of filtration of the water, I& is the permeability in the two-phase zone, and pWois the density of the water when there is no excess pressure. The permeability of a porous medium is related to the humidity v by the following relation that is an approximation from experimental data [3]:
Moreover, it is assumed that the thermal capacity and the thermal conductivity in a two-phase zone obey the mixing rule: plllC~-(i-m)p,SC.~+mvp~~+m(l-v)p~i, L-(l-m)h.,+mvb,+mji-v)t
ppt C, and & are the density, heat capacity and therma conductivity of the surface of the ground respectively. The boundary conditions are integral analogues of the equations of conservation of energy and mass of the water, recorded for the boundaries of the zones taking into account that at the boundary between the liquid phase zone and the two-phase zone L’s(t) the humidity Y returns to I and at the boundary between the solid phase and the two-phase zone Z,(t) it takes the required value Y.:
Here
r-Xl(t):
v-v.,
T_-T+-T.-To-tjP.,
The initial conditions take the form T@, 0) -To,
P(z, O)=O,
T&f(O)-Th
x, (0)-X*(O) -0.
(11)
If in the boundary and initial conditions (10) and (11) the quantities To and 7@are constant then the problem has a self-similar solution of the ordinary type: T-T(E),
P-WE), x, (1)+‘I,
v-v(E), XI(i) -IW.
+rr-‘h,
(W (12b)
This reduces to solving a boundary value problem for the system of ordinary differential equations that have been set out clearly in (7) for the solid phase domain (KKy) and in (8) for the liquid phase domain f&C%=) for Tr=T,-#.
94
The constants y, /3, P’ and T. have to be found from the conjugation condition of the solution through the two-phase state zone in which after substituting (12) into (9) we have h,T”+T’v’(h,-hr)nr + $
IT’p,C,E+mTv’(pC,-p,C~)E+p,m~v’El-
0,
(13a)
At the boundary [=y the conditions I-T.-T,-qP.,
P’=
\-
A,T_‘-A,(v.)T+’
- +
mt~p,r.y,
are satisfied and at the boundary [=& the conditions T-T’-T,-qP*,
!!=i,
T-‘-T+‘,
P-‘-P+‘.
(15)
hold. The non-linear boundary-value problem (13x15) is solved numerically. In Fig. 3 the results of the calculation of the distribution T(the solid line), P (the dashed line) and Y(the dot-dash line) are given for the following values of the soii parameters [4]: Z&274%, ?%?63%, a,=0.7gX10dmz/s., u,~l.14X10dm2/s, ~=0.86X10“mz/s, nH.2, k,=8.SX10‘“mz, pS=2000kg/m~, C,=0.92kJ/kg°K, A,l=2W/m°K.
\f ._.-----. 6 Flu2v
\
.A-\
g
1 ‘1, ‘.
1s
6 06 4
.4
FIG. 3
The results show that the frontal model usuahy used for engineering caIcuiations can lead to an inadequate description of the freezing process of a material of low permeability. The realization of some or other freexing state depends qually on the value of the Permeability and on the values of the initial and boundary temperatures. These Parameters can be represented in nondimensional form as follows:
FIG. 4 In the plane of the Parameters (VO,& (Fig. 4) the position of the boundaries of the domain of realization of the freezing regimes for diierent values of $ is shown. The domains lying below the curves correspond to solutions with a phase transition front and those above the curves correspond to values of the parameters for which a two-phase zone is created. A decrease in $ leads, of course, to a widening of the domain of existence of the two-phase zone. Trandared by S. R.
I
95
REFERENCES 1.
TIKHONOV A. N. and SAMARSKII
A. A., 77~ Equations o~~azke~ticu~Phy~ic~
(Uravnenie ~temati~b~koi
fiziki), Nauka, Moscow,
1966.
2. 3.
4. 5. 6. 7.
COLLINS R, Flow ofF&d Through Porous Muteria~s, Mir, Moscow, 1964. Thermophysical Research on the Cry&e Zone of Siberia (Teplotizicheskic issledovaniya kriolitozony Sibiri, Nauka, Novosibirsk,1983. KUKHLING KH. Han&ok ofphysics, Mir, Moscow, 1985. IVANTSOV G. P., “Diffusion” supercooling in the crystallization of a binary alloy, Dokl. Akud Nauk SSSR, 81,2, 179-181.1951. BORISOV V. T.. The crystallization of a binary alloy with conservation of stability, Dokl. Akarl Nauk SSSR. l&3,583-586, 1961. ABDONIN N. A., A ~athe~ticuI Description of the Processes of Crysta#ization (Matematiche~oe opisanic protsessov kristaiiiaatsii), Zinatne, Riga. 1980.
(NW-5553/86 glO.~.~ 0 1988 Pergamon JournalsLtd.
U.S.S.R. Compztt. Maths. Moth. Phys., Vol. 26, No. 6, pp. 95-98, 1986. Printed in Great Britain
EXACT ESTIMATES OF THE ERROR GRADIENT FOR PROJECTION-DIFFERENCE FOR PARABOLIC EQUATIONS IN AN ARBITRARY DOMAIN*
SCHEMES
I. D. TURETAEV Estimates are given for the rate of convergence of some two-layer economic projection~ifferen~ schemes @.d.s.) and Crank-Nicholson p.d.s. for solving parabolic equations in an arbitrary domain. 1. An optimal estimate of order O(r + AZ)is obtained in [1] in the norm of h(Q) for the error of a family of economic p.d.s. for solving parabolic equations in an arbitrary domain fl c R*. In the present paper we prove for these p.d.s. orderwise exact estimates of the rate of convergence O(s + A) in the norm II,II(*J-IID~.R~+IID.u*. _. In addition, we study the Crank-Nicholson p,d.s. in an arbitrary domain fl c I?“, n Z 1. For the error of this scheme we prove exact estimates order order O(# + ]k]*) and O(S + ]hl) in the norms of f&Q) and ll.@r’respectively. For the case of constant coefftcients, when ff is a parallelelpi~, the optimai estimate of the error of the Crank-Nicholson scheme in the norm of L,(Q) was obtained in [2]. The rate of convergence of the p.d.s. whenfsatisfies some smoothness condition was studied in [3,4] (see also the references there). We solve the first initial value problem Dlu+Lu-f
in
Q-$2X(0, T),
ulho-uo(+),
u/r-o,
I%%aX(O. T).
(1)
where ~~-D~(Q~~(z)D~~) -+I u,r=Q (summation from f to n is understood with respect to repeated subscripts i, j), and fl is an arbitrary bounded domain of R: n L 1, with dfIcC. We assume that ~l~i*catrE,E,~r~lll* for any f = (51, * . . . 5.) = R, &%I (y&r, 1.1is the eikonal of the norm in R”, and o,~-o~(, an=W-‘(Q), o=&.(R), 6 I -4, 2,. . , n. Under these conditions we have for the solution of problem (1), (2) the estimate llD,Pull,+liD*Dtull~+llLII),nli~, ,
Here, Drp is the gradient of 9, D29 is a matrix with
D,Dw, ndil~-HD~f~l*+iif~n, +ILc& elements (here, are the norms in L*(Q) and f,Dcf=&(Q), fo, .a,&=U+(Q); Il.ilp. 11.11, @*(n) respectively), and K is a constant, independent of the functions in the inequality (and of r, A, if the data of the problem depend on them). We assume below in this section that n = 2. We construct a mesh domain (see [ 1, 51). Let R? be a rectangular mesh in R2 with steps A,, !tz with respect to the variables x,, x, respectively. Let &, be the domain which is the union of partial rectangles of Rh2which lie in fi and are spaced 0 (max (h2, !I,)) from the boundary XL Further, the domain a\& is divided into triangular elements in the usual way. Then, the required mesh domain f& is defined as the union of all triangles and f&. We assume that the division of fl into finite elements is regular. With this way of constructing fh, the boundary dfl is approximated by the boundary dfl, of the domain f& with order O(K), where A is the maximum of the sides of the elements. Let s” be the space of functions, continuous in fi, linear in triangles, and bilinear in the rectangular elements, and vanishing in ii/n,; S, is the space of functions, continuous in [O, r], linear in each interval fr(k-tf. rk], k=i, 2.. . . ,.4&TM-T. We put S = s” X S. We consider the following family of schemes for the approximate soiution of problem (1): (yr”‘, cp)+‘WDy?“, +a‘acym,
‘Zh. vychisl. Mat. mat. Fir., 26,
11,1748-1751,
1986.
Dcp)+O+(&Dxytm,
DiDtcp)p~
cp)==(~fl”. CPL m-o, 2, . . . ,Y-i,