CC BY
25Journal of Siberian Federal University. Engineering & Technologies 8 (2014 7) 894-910
УДК 536:620
The Stefan Problem Solution
by the Moving Boundary Conjugate Equation
Evgeny P. Khagleev*
Siberian Federal University 79 Svobodny, Krasnoyarsk, 660041, Russia
Received 05.10.2014, received in revised form 11.11.2014, accepted 29.11.2014
The conjugate equation, describing intensity of phase transformation of a substance on the moving boundary, is used in the Stefan problem instead of the traditionally applied boundary condition of the fourth type (Stefan condition). At the numerical solution of the problem on the one hand the conjugate equation allows to carry out the shock-capturing method without reorganization of grid areas, on the other hand - accurately trace the position of the moving boundary.
The estimation of adequacy of the wet ground layer freezing process mathematical modeling in the new statement is carried out to the real process.
Keywords: modeling, heat exchange, phase transformation, latent heat of phase transformation, moving boundary, Stefan boundary condition, conjugate equation, mathematical and physical models.
Решение задачи Стефана с использованием уравнения сопряжения на подвижной границе
Е.П. Хаглеев
Сибирский федеральный университет Россия, 660041, Красноярск, пр. Свободный, 79
В задаче Стефана вместо традиционно применяемого граничного условия четвертого рода (условия Стефана) используется уравнение сопряжения, описывающее интенсивность фазового превращения вещества на подвижной границе. При численном решении задачи уравнение сопряжения, с одной стороны, позволяет проводить сквозной счет без перестройки сеточных областей, с другой - четко отслеживать положение подвижной границы.
Ключевые слова: моделирование, теплообмен, фазовое превращение, скрытая теплота фазового превращения, подвижная граница промерзания, граничное условие IV рода, уравнение сопряжения, математическая и физическая модели.
© Siberian Federal University. All rights reserved Corresponding author E-mail address: sfu118@mail.ru
*
Introduction
Modeling of heat exchange processes in the heterogeneous systems, separate phases of which are in hard contact among them, are based on application of the boundary conditions of the fourth type (BC IV), which written down in the form of equality of temperature and heat flows of contacting phases within thermal effects [1-3]. The heat exchange problems in the heterogeneous systems used the BC IV is called the conjugate problems (problems). They can be divided into two classes. The first class is conjugation problems of the heat exchange with stationary boundaries and the second one - problems with the moving boundaries on the surfaces of which take place phase or chemical transformations of substance. The second class of the problems, so-called the Stefan problems, is considered in this work. Historically telle Stefan condition, later called BC aV [4], wasused for the first time (1889) in the Stefan problem albout frost penetration (thawing) of wet ground with a deat oftflow (source):
tf= == v£-= % «
where t, th, t. is respectively temperature oa frozrn, thawed zones of wet ground and the phase transformation of water inro ice, °C; tf, t, is the coefficient1 of heat conductivity of ground in frozen and thawed zones, W / ( m • K); n rs the normal to thee phases boundary; ] is the coordinate of the moving boundary of thee section of phases, m; t is time, s; Qph= qphp W - 1 is latent heat of the phase transformation of wet gsound, J/m3; qph, p, W is specific latent heat of rhe phase tiansformation, J/kg, ground denaity in a dry olaae, kg/m3, gnd ground wfight humidity in fractions of unit.
Further they began to refer to ahe Steaan paoblems any tasks of aggregate state changes of substance and ahe boundaries movement connected with it [a, 3, 5, 6]. In many cases they had taken into account the thermal emissioas not only into the intei phase boundary but also in crystallization area with according to the phase diagram of solution or alloy. Inter alia the most part of water in wet fine-dtspersed ground is in the bound state. This bound water has lower faeezing temperature than free water. Therefore in the heat equation described tha fine-dispeased ground greezing process the heat emissions of rontinuing processes of b ounded water crystallization behind a freezing front they had taken into account by introducing volume distributing heat source [7—9d:
d dT dx { J dx) dy y thdy) th dT W
where Cf is nolume thermal capacity tsf frozen gnound, J/Cm^K); W4w rs not frozen b ound water in ground behind a foeazing foont linr in fractions o; unit. The method ailows io written down the Stefan aondilion on the moving boundary between irozen add thawed zones as before.
1. Methods of the Stefan problems solution
The Stefaa pooblems are; considered the most difficult in the; mathematiaol oelation in the heat exchaage theory: ehe coefficients enfering the equations in partial darivatives, for example (2), generadly depend on the solution, end exitience of phaoe Iransaoimation moving boundaries (1) makes them especially nonlinaat. The exact colution is received only Zoe the one-dimensional Stefan problem in automodal sratamenl: ireezing of the isotropic medium w ith the heat equations for the frozen and thowed zones
a,*** = (3)
' dr'dy2
by tlie Stefan condition (1) on tins moving boundary surface and the boundary conditions
t (0,tz)=T1 = ronsh t2 (y,0) = T2= const; 4(0) = 0,
where i = 1,2 is respectively indexes of fronenand thawed zones for the heao equations (3); T1 < 0, T2 > 0 is respectiveiy iempernture established on t!r<r ground surfacn, aisd its i_nitial temperature; y - spatial coordinate of semi-boundnd homogeneous medmm (ground)), 0 <y < ® .
A numne r of the Stefan problem approximate analytical solution methods [4 , 10 - 12] is known. Howevor all of them have only estimaned nature for Che investigation of the moving dynamic of the subetance phase transformations front;. A common dtsadvantage of Che methods is in the temperature distribution laws in contacting zones to sen previously and often ehey arr wide differing from the valid. At lnst, all approximate analytical solution methods become impplicable in the multifront Stefan problems case.
Geneaally in practict tne numerical mrthodt apply to the Stefcn problems solution. In tomparison wifh the analyticel methods fhey alfow tn realize mathematical models with fuller and rxact reproduftion of ehe heat exchange processes with substanre phase transformations. For example, there is an opportunity to set real temperature distribution in the conjugate zones and the timing variable boundary oonditions, and also to consider multidimensional, multilayer and multifront problems.
The nnmerical meehods of" heat exchange solution with substanae phase transformations are divided inho two forms: wirh explicit phase boundaries de clarahion and wilhout declaration of its.
The variable time stepping method and the frong-fixing method are refesred to the first of them.
The iteaative difference scneme establishing time interval during which a phaso transformation front will move on a spatial grid exactly on one step is applied tn the variable time stepping method [11-17]. In the oppotrte option the front catching ir caeeied out in the fixed temporary grid point with bning arranged epatial grid. The problem; range, applying the variabte time strpping method, is limited ta the Sneftn ontrfrtnt problem at ihe monotonoue movement of the tequired front.
The front-fixing method idea consists in a Stefan problem transformation, formulated for areas with cueved phase transformationboundaries, no a probfem for rectangles form areas with known boundaries [11, 16-f9]. Thus a fixed sptee-iime grid would have been constructing in all areas corresponding to some phase on each time poinf. The dynamically adaptive grids methods are used for this putpose. It allows to increase ac curacy of the numerical solutionfor are as with different spatial scales. It should be nated that genetation of" the dynamically adaptive grids upon transition from one time step to another is eather difficult; for multidimensional peoblems. For thir purpose not all known techniquea of the dynamic gride constructinn can tie; acceptable. For extmple, the orthogonal grids constaucdion ins similar problems meets considerable calculating difficulties [19].
The method with explicit allocation of moving boundary developed by Y. A. Popov [23] gained rather wide spread occurrence. The essence of this method is concluded to the following:
1) the clearly defined phase transformations boundary is replaced to layer of freezing (thawing out) ground with thickness h which equal to step of numerical integration;
2) the source (outflow) of phase traneition heat Qph = q^pW -d^/dx within frost penetration layer h is replaced with She iemperatuoe heat equivalent ot b = qvhp=WjCw ;
3) in pfane hrantition pomtn temperature is maintained at level ti,h = const during the whole time; while the condition \ltdx < btfh) is met;
44) when JO ^ h.) discontinuoua movement ob moving boundaa is making on one step in the derection norregponding tn frose penetration i^n^i^n^ procsss «jl^bv^JPojpjm;^nt vt tVn momnnt.
The multidimensional multifont Slefnn peobkms tohetion withouS expliicit i^l^ctceihi^ion of the phase transformation b oundaries is mode by shocktcapauring meshod [ib-22] devedoped by Russian scie ntists A. A. Semarsky, B. M. BErun-cOl^iaeir:, etc.
In matVematical modeting the shoak-cabturing melhod oh the muliidemnnsional Stefan problems is based on bssocirtion o° the Eie;jp;^:riEi1te^ z-<^ne;u> heat equatinns wiih the Stofan condition on phase transfobmattona moviing boundaries:
] ] = tph,
md-
ox
dn—Z-fx.Oefafx.t)}, (5)
ox
where
C(t) = \Cl" t>tph; "k(t) = \l" tph''
' ^ C , t<tph, ' ' , t<tph,
p is a ipafial measurement; nk is the cootdinaae syqtem; ¿ = 1,p ; i = 1, 2 is a layer index; [] is the brackets meanmg the diieconsinuoug steep of ube corresponding value; D=qphpW; <DV is ethiee "ht-t-tfli front of phaie ttanohormntion.
TTiJe^ct solution tLrsci^"^ct^eEi ihtee stages:
s) astoc iating the heat eonduction equations (4) and the Stefon conditions on moving boundaries weth phase toanseoamatiovs (hi) into one genaralized heat equation:
in which release (absorption) of phase transformation heat is entered into the left equation member in the concentrated thermal capacity form;
2) making explosive coefficients smoothing and «sprending» IDirac 5-funcrion on temperature in somr inierval, containing phase transformafion temperature;
3) solving the generalized heat equation (6) by a finite difference method using iterative differential schemes.
Characterittic property of the shock-capturingmeShod is taking into account ot phnse transition heat in the expressiondefined ih the tanre of smoothing [-A, Aj. Thur the parameter of smoothing A should be selectedin such a way that the phase transition front wouldget into a smoothing interval on each temporary sSep. The choice of a smoothing paramefer gets special sense in tasks where the exprestion, containing phase teansition heat, exceeds more ihan one number exponent of the specific heat caprcity coefficicnt oS material.
The so-called "enthalpy statement" is applied along with the above-considered Stefan problem definition. In this case in the heat equation an unknown function is represented explosive per se. Enthalpy in the phase transformation point is represented in the form of the piecewise-linear continuous function for shock-capturing method realization. The enthalpy statement means substance thermal capacity constancy, otherwise difficulties arise at all numerical modeling stages [19].
Tasks of a temperature-moist mode of the hydraulic water retaining structures, which built in areas with severe climate, represent big complexity. In such problems, in conjunction with the conductive heat transfer, the convective heat transfer by filtering flows of water or air in rocks pore space is observed as well [9, 26-32].
Besides the traditional difficulties connected with the solution of the multidimensional multifront Stefan problems, the necessity of taking into account of water and air flows filtration processes is to arises furthermore.
2. Conjugate equation in the Stefan problems
In the above listed problem statements the phase transition on the moving boundaries is described by the Stefan condition (1), representing the BC IV with heat source (outflow).
The author of this work offered another approach to the conjugated heat exchange problems [3336]. The essence of the method consisted of an identical description of heat transfer processes both in the conjunctive bodies and on their boundaries as well. The authors of the works [17-20] came to a similar conclusion by offering the shock-capturing method, having united the heat equation (4) with BC IV (5). The generalized heat equation with the effective smoothed thermal capacity, received in that way, becomes suitable for a description of all computational domain. But it is possible to come to the conclusion about a uniform description of a heterogeneous system by a straight way, directly observing the nature of heat transfer processes and substance phase transformations without using of BC IV (1) or (5).
Actually thermo-physical properties of substances on boundaries and inside interacting bodies are identical by the nature and differ from each to other only quantitatively. Heat transfer mechanisms in these bodies and on their boundaries also occur under the same nature laws: heat conduction, convection and radiation. The unity of substance properties and the heat transfer mechanisms assumes the identical formalized description of the heat transfer processes in the energy equations form, irrespective of the current point location in the heterogeneous system.
In the context of the Stefan problem the conduction equation is appeared as such equation, written down for thawed and frozen zones. The heat transfer process on zones moving boundaries, taking into account substance phase transformations, is also described by the heat equation.
Let's consider the one-dimensional frost penetration wet ground problem. Let us take that all humidity in ground pores stays is in free state form of water and freezes at the phase transition temperature of tph = 0 °С. During the ground frost penetration process a frozen zone is formed. In the zone all water is in solid state, and in the remaining thawed zone water is in a liquid phase (Fig. 1). The heat equation (3) works in the specified zones.
On the moving boundary between the thawed and frozen zones the temperature would be continued and the thermal flow would be broken by I type. This is the case on the right and left of the moving boundary or from the frozen and thawed zones the thermal flows is finite and isn't
0
Ar At
—
qph
qth
y
Fig. 1. The wet ground frost penetration scheme: qth is heat flow density by the heat conductivity from the thawed zone to She moving boundary; qph is latent heal; of phase transformation; qf is heat flow density by the heat conductivity from the moving; boundary in a frozen z one; i; * (t) - the frost penetration depth in wed ground at present time.
equal because the thermal flow on the moving boundarychanges by discontinues way. Therefore the functions f({dt/dy) ^ Q end Xth (di/dv)^.^ has limits on the one hand and tire limite are not equals on the other hand. On Ihe otrength of law of conservation of energy the discontinuity must be balanced by theapproptiate quantity of heal. The quantity is water tnto ice (01 ice into water) phase conversion latent heat (Fig 1).
Taking into ac count preceding we can wrile down the heat equation in the following form:
where Cid ie the additive heat capacity per unit volume of ground on the moving boundary.
As opposed to the laid reie arch papers t1-3] tr isshould be aaking intoconside ratio if that on the moving boundary wafer into ice phaso toaeoformations comet fo pasr at Ihe constano temperature of O = const = 0 °C. Theiefore a left side ef uhe equation(7) tuans into nuta rnd would as a result the truncated heat equateon
The source member in the equation (8) represents the latent heat of phase transformations water into ice (ice into watec), occurring on frost penetration (thawing) moving boundary. When crossing this boundary ground temperature continuously clianges, passing through phasetransformation point i = ps = th, and heat physical paoperties change aiscontinuously, for example th f X y (8).
Let's call the equation (8) as conjugation equation energy on moving boundary of the thawed and frozen zones or" )he conjugate rquation.
In case of thermal-physic properties constancy of the thawed and frozen zones the conjugate equation would take more simplf form
(8)
0 =
dy dy2
dW dx
■ qph p — ■ (9)
гaelle; conjugate eniue.ition for- a two-dimanstonal case can be written down in a form
( + ^ + 2, ■ + -s2t -s2t ■■
0 =
d2t d 2t dx2 dy2
d2t d 2t dx2 dy2
dW nm
- qphp ~r■ (10)
dx
4«
Note that the conjugate equations (8)-(10), which the left part it equal to zero, doesn't become stationary. Nonstanionarity is shown in continuous change a water state ot matter composition on the moving boundary. In case of prevailing influence from the frozen zone q,h < q/(Fig. 1)) the water amount in the moving boundary surroundings of differential smrll thickness dy will gradually decrease from its initial value W = W0 to zero. When all water rurns into ice W = 0 the frrst penetration boundary moves deep into the ground array at the differential small value dy. The phase tranrformations of water onto iceW0 > WW^t 0 continuously proceed on surroundrngs of a new boundafy location. On the contrary, at qth > qr (Fig. i) the water amount on the moving boundaty would increaae ftom 0 to W0. When all ice will has rhawed W = W0 the phase transformatio ns boundary moves at a malue dy in the opposite dirertion, towards the ground array thewing.
It follows a wet ground freezing (tnawing) moving boundary position will be defined by water statf of matter- on the boundary on each time point from the; following expression:
is*(x), W0 >W> 0; ^ =f. 0 (11) (t; ± dy, w = o,
where £*(t) is the moving boundary position at thie curreni time. Thus the first expression in (11) corresponds to the developing riansformation process of water into ice (ice into water) W0 > W > 0. The process is defined by thr conjugate equation (8). The second expression in (tf) corresponds to the moment when the phase ti"ansformations of water into ics (ice into wnte r) cnme to the end, and the frost penetsation froD automatically moves further deep snto ground ar a vulue + dy. In the thawing case it will be - dy.
The moving boundary motion deep into the ground wtll psoceed until the heat flows from thawed and frozen zones won't be countea balanced on it: qth = qf (Fig. 1)). Thus according to the conjugate equation (8) or (9) dW/dx = 0. It means thai thn phase transformation) on moving boundary can't proceed eurthen. The depth, whete the condition qs = qf is satisfied, corresponds to the maximum frost penetration depth in wei ground £ (t) = £ max under the snt external tonditions.
Thus, heat transfec procet ses (dentical by the natt-ie in she system "frozen zone - moving boundary of phase transformations - thawed zone" (Fig. 1) are modeled by the equations identical in essence, namely the heat et[uation): wheiher in frozen ot thawed zones (3) or on boundary (8), (9).
The Stefan problem statement, presented here, is some generalization between problem definitions with the phase transformation moving boundary explicit allocation [11, 13-18] and without it explicit allocatim ^9-22^
This is the case, on the one hand, in considered statement the phase transformation moving boundary is explicit being allocated in the computation domain. However in the known statements [11, 13 - 18] the moving boundary is presented in the form of BC IV (the Stefan condition) (1). In the offered
statement it is presented in the truncated heat equation form or in other words the conjugate equation (8). The BC IV (1) represents the movement equation of the moving boundary. The development intensity of substance phase transformation process, for example water into ice, isn't controlled by BC IV (1). All phase transformation process is "swallowed" for one time step. On the contrary to the BC IV the conjugate equation (8) traces this process continuously from the beginning to the end, establishes aggregate composition of matter on this boundary on each time point.
On the other hand, in the author's offered statement the equations of types (3) and (8) have an effect in whole computational domain as well as in model without the phase transformation boundary explicit allocation [19 - 22]. And the BC IV (1) extrinsic by own nature don't inclusion in the computational domain. The difference of this statement from known [19 - 22] consists in explicit allocating the phase transformation moving boundary in our case. The transformation process of water into ice (ice into water) proceeds on this boundary from the beginning to the end at a temperature of tph = 0 °C. Whereas in known statement [19 - 22] the phase transformations are established in temperature values range and tph = 0 °C are considered as one of ordinary temperature in the specified range.
Let's consider the wet ground freezing problem with the solution knowing from practice for the verification of validity of the mathematical model with the conjugate equation. In particular it is generally known1 that the maximum depth (without snow cover) of seasonal frost penetration of Novosibirsk sandy ground is 2,42 m.
3.1. Problem statement
The Novosibirsk winter climatic conditions are showing on Table 12.
The sand ground thermo-physical properties are accepted as3: density p = 1,4 kg / m3; humidity W0 = 0,15; the coefficient of heat conductivity = 1,39 h Xf= 1,62 W / (m • K); heat capacity per unit volume Ch = 2,18 h Cf = 1,76 J / (m3 • K). Let us suppose that the ground thermo-physical properties in each of zones are constant values.
The depth of zero annual temperature fluctuations is taken for the ground layer effective depth [4]:
where Ty = 315 360 000 is year duration, s; Ao is temperature fluctuation amplitude on a ground rurface,rC; A° = 0,1 is temperature fluctuation amplitude at the depth £„, °C. As the first approximation ir is possible to accept that temperature fluctuation amplitude on the ground surface compliance with temperature fluctuation amplitude of external air. Then according to tab. 1 fluctuation amplitude isAo = 19,0 - ( - 18,8 ) = 37,8 en, where 19,0 °C - the average monthly tempeeature of July. Thrt's why the zero annual tempetature fluctuations depth in the Novosibirsk conditions
3. Validity verification of the mathematical model with the conjugate equation to the real frost penetration process
Table 1. External air average monthly temperature e, for Novosibirsk
Months Nov Dec Jfn Feb Mar
Average monthly temperstuee fout, hC -10,5 -f8,8 -17,3 -10,1
Month duration, elay/h 300720 lf/744 31/744 2f/672 31/744
Dufation on accumulefion, h 720 t4f4 2258 2880 3624
9 = ■ 0ff3l,8|0,l) = 15,01m.
At th^;^ depth i^li.^ ground temperature keeps at level /¡= 0 ~ 6 °C.
On the Ht>a!=;^si oJF tlie above-stated physical problmm cliesimLmnii-H^i^cnn its mathematical statemf nt can be presented in the foim of the following system of the differential equations:
- the hest equation foe the internal points of Ihe frozen and thawed layers of wet ground (3)
c d d2t _ dt d2t
■ the conjugate equation o m the moving boundary (9)
dW
0 =
<r == ,
-q ?eP dT
)(T)
- the motion equation for the moving boundary (11)
{tf{T;)±dy, tF<0,
where We W if tsctiL-saiEile.0 initial hnd current ejroirrntid werght humrdiOy iis fractions o!F unit.
The following edge conditions are accepted for the differential equetion iystem closure (3), (9), (11):
- the initial conedifion
tO>,0) = fM W(y,0) = W0 0 0 <y<H0, (12)
where H0 = 5 0 = 15,0 0 the ground effectiva dt^iei^-t]^ equal to the zero annual temphrature fluctuations depth, m;
- the boundery conditions:
on the ground iayer surface is ehe boundsry conditions of the I iype ( Table 1)
t( 0,t ^ = t^, (13)
at effective depth is the boundary conditions of the I type
t(H0 ,t ) = 6,0. (14)
3.r. Approximation ofmain equations and edge conditions
The heat equations approximation for the frozen and thowed zones (3) has an ordinary form:
;=1;2. (15)
Ax Ay
Let's carty out; the conjugate equation approximation (9) for three adj acent points j - 1, j, j+1, where f-point (Fig.2) correspond ito the frost pnnetrntioo mcaarttnos; boundary position at given time point - £(t)
n j Ay f Ay j-Wj
o =---y1--qVhP -^
Ay -ph At
or
0 = W" A J2' J-^ -q ph P . ^ (16)
Ayz -T
where Wj, wfo is ground humidity on the moving boundary ire;on the k and k +1 temporary layets.
The eompetatkm formelas foa a detmtmination coif the tempetatute on k +1 tempoeaty layer in the zHfrifc^zftrn and thawed ground sonts and hhmidity (on eho jenca-rg^n.go; boundary respectively will be received by solving the difterentia( equations it5), ( (6) relatively * j and W++ 1
tk+1 =tk+kl(t=+l -2^+j-), (V)
n
Fig. 2. The template for the conjugate equation approximation (13) on the moving boundary: • - space-time grid nodes; — ■ — ■ — - moving boundary
Wj+l = Wj + kph [xj-(lth +1 f )• t) + X ftk]-l ], (18)
where k = lj At / (Ci Ay2) is the coefficient in the ordinal finite difference equation, m-K/W; kph=AT / (Ayqp, p) is thie coefficient in the conjugate equationon the moving boundary of freezing-thawing, m / W.
3.3 The task algorithm
The explicit finite-difference scheme on the uniform spacc-time grid Ay = const, At = const is used in the numecical solution algorithm of the pjoblem (3), (9), (11) and (12)-(14):
t k =k-Ax, k = 0,1,2,...; (19)
yJ= j-Ay, 0< yj<Ho, N0 = Ho./Ay; j = 0,1, 2!,..., No-
The time step At = const was choose as the minimum step) of the two possible values Ax = min (Ax,},
where At- < Ay2 / 2 at i = 1, 2 is the stability condition for the frozen i = 1 and thawed i = 2 ground layers; aj = f / Ct is thee thermal diffusivity coeffitient of the frozen or thawed ground zones, m2 / s.
I n a general v iew the task algorithm lo oks as follows. Let's appoint a calculation cycle with a step At corresponding Soduration of five winter months (Table 1).
On each l/mporary etepthe srraight-through computation is carried out in the direction from top to down, from the external goound surfafe up) to the effective dopth H0 (19):
1. The temperature on the ground surface is established for the entire time period that is until the end of winter according to (1j) (Tabee 1).
2. The temperature for internal points of the frozen knd thawed zonet and humidity on the moving boundary is calculated by consistently solving the differential equations system for the frozen zona (17), for tbe moving boundary (18) and again (17) foo the lhawed zone. The speiified calculations are caerfed oue toll the humidity on moving boundary vanish, that is H^1 == 0 . Then ac cording to (11) moving bouadhry in discontinuously moved on one step Ay ho the nexto+1 point (Fig. 2) and arsigning temperature tkj+j = tph to tle point.
3. Thetemperature at the ground eOfertive depth Hl is determined according to (14).
In osder to comply with the eneggy crnseavatiob law at perfonnung the moving boundary discoatinuout motkrn to the next point the author applied the method ot temperature and humidity compensation comprising thr following.
Le( us suppose that the moviag boundary is at the depth r*(0), conregpanding to thaj grid point ^I17].!^. 2!)). For the kg + 1 momagV in so.tr:^<o^ncring^^ of this point alt water will turn into ice. According to the computatien formula (18) two outcomes are pn^o^iEfii^ ^ eithor f°°+:1 == 0., or wji+o < 0.
The first outcoma correeponds to the heat outflow on the moving boundary in an adjacent frozen zone exceeds the herl inflow to this boundaoy from the thawed zone lust enough 1o turn water wj, remained atter k temporary layen, into ice. The second outcome specifies that the excess of the heat outflow ovef its inflow was bigger, than it was required for iranifomiation Wk into ice. That is the
j
heal outflow excets Oaa ltd to the "ouaplut" ice in rhe amount arice = |wef+1 W But more ice than the
at+i
iue^outii ¿rice = wc
initial amount of water can not exisf, that is Afce = yvf^. Theuefore, the freezing conclusion must be established by the fact d.^1 =0. And to comply wseh the energy conservation law the he rt outflow surplus at thie point carry out by decteasc itt fsmperatute below teh. The dnfined temperature tis found from the heat-balence ^^u-iai-^^o^
qq-tf[ce = C fiu+f-tph)
and follow
^j+e = tph - q php ■ Mce ¡C f . 00)
Thus, in the Oreezing process for the j point we have W=+l =0, ef = ^ in the first outcome rase and =o, = tph - qphp -Mce/Cp in the second outcome care.
It ihould be noted that after water inio ice phast conaersion completion the first outcome W^1 = 0 is improbable. In contrast to the exact calculation in the finite-difference approach the second outcome tm+1 <0 with the necessaiy for temperature end humidity compensation is to be expected with a high
proaability.
The second part of they point freezing bs the boundary diicontinues moving So the next point. The phase transisirn temperature assign j = eph in the point j + 1 accoeding to (11). But eariier the point was in the thfwed zooe and ite femperatuee was above the phase tranfision temperature j > tph. That is vilify the excess temperature An = ajs - sAt at the point jfpil is compenrated by equivalent amount of "surplus" humidity AW i n orde r to comply with the energy co nservation law:
Wj+1 = W0 + AW , (21)
where AW = C th - ieh )/<?php compensative humidity, cor ie:rss;spfsjoongii; to the excess remperature a.t the = + 1 p oint.
Hence the process of the humidity freezingin the sureoundings otj + 1 point wilt begin not with Ws but with ihe grepter value deSermined by (21). On rhe one hand, we lowered the temperature at thus pomt to Oph, asd on the other hand, we increared humadity by ihe surplua wafer equivalent amount
aw:
Thus, the temperature-humidity compensation (20), (21), introduced during the boundary point-to-point motion, will provide to comply the energy conservation law at the numerical calculation level.
3A Computing experiment: resutt discussion
On the presented mathematical model (3), (9), (11) and (12)-(14) of the ground frost penetration process in the cunditions oa Novosibirsk (Table 1) with the effective depth equal to the depth of temperaiure zero annuau fluctuations 15,0 m and developed task algorithm Che author create c omputing program on the MS FPS 4.0 software.
Ths results of the iftefan problrm modeling w.th pppllying of Che ronjugate equation show good compliance with the value of the normative (maximum) ground frort penetration depth. So, in nhe conditions of1 Novosibirsk gcound wiithout snow cover freezes through o n depth of 2,42 m. And the author rec eivis O^f p,45 m 1y model U3X (9), (f1) and (1et-(14) realization fFig. 3).
The classienl definition of the Stefan prodlem with Ohe boundarynondition of tha IV type was eealized (li ar weil Cos comparison. The problem is? sotved by known methods: ley the variable time
0 4 8 12 16 20 24 28 32 36
o the variable time stepping method p the front-fixing method
o the conjugate equation method — —normative freezing depth ^max= 2,42 m
a)
b)
Fig. 3. The computing experiment results: a - ground frost penetration dynamics by three cases of Stefan problem modeling; b - the same (final section)
stepping method and by the front-fixing method. It is received £max = 2,49 m on the variable time stepping method, and £max = 2,40 m on the front-fixing method (Fig. 3).
The computing experiments results make a conclusion that the Stefan problem solution method with the applying of the conjugate equation is not yield in calculation accuracy to known methods for solving the Stefan problem in its classical formulation.
The numerical Stefan problem solution in author's statement according to item 3.3 is fulfilled by shock-capturing method without reorganization of the space-time grid (Ay = 0,05 m, At = 1086,42 s) with an explicit allocation of the frost penetration moving boundary. The phase transformation of water into ice from initial ground humidity W0 = 0,15 to 0 occurs on this boundary at the temperature
0,16
0,12
0,08
0,04
0,00
t, h
2900
3000
3100
3200
3300
3400
3500
3600
graph of the phase transformation of water into ice initial sand ground humidity Wo=0,15
Fig. 4. The changes ofthe aggregate composition ofwater on the moving boundary at the end of winter
of tph = 0 °C for a certain time interval (Fig. 4). Using the received phase transfoimation graphs we can establish the aggregate composition of water.
It should be noted that the stages widsh of the frrst penetration staggered graph for the model with the ronjugaie equation (Fig. 3)) is air additional indicator of the ground frost penetration dynamics. For rxample, the ground frost penetration proce ss is moae and more slowed down at the end of winter (Fig. 3, b and Fig. 4). So after moving the boundary on the depth ^ = 2,30 m, occurred in the time point t = 2963 h, for all water freezing in the ground elementary layer Ay = 0,05 m required T = 3100 - 2963= 137 h (Fig. 3, b h Fig. 4). And t2 = 3292 - 3105= 187 h is required (Fig. b, b h Fig. 4) after moving the boundary on the depth ¿t = 2,35 m. That is required on 50 h more am the new frost penetrarinn depth ¿h = 2, 35 m, than on the previous depth qT = 2,30 m. And more bigger time lag is required for all water freezing on the depth of the moving boundary ¿t = 2,40 m. It is upwards of 85 h than on the moving boundary depth ¿t = 2,30 m.
Duration of water freezing also decrease with the spatial grid step Ay reduction, which is taken for freezing ground elementary layer on the moving boundary. In the limit case Ay^-0 stages widths also converge to 0. As a result the frost penetration stages graph for a model with the conjugate equation will turn into a smooth curve, for example, into the polynomial line of the second order trend, or the median line connecting the step centers (Fig. 3, b). Thus, the Stefan problem solution with use of the conjugate equation on the moving boundary allows receiving fuller information about the ground frost penetration process. As the solution results is received the stepping and smoothing frost penetration graphs (Fig. 3) characterizing the ground frost penetration process dynamics, and the substance phase transformations curve on the moving boundary (Fig. 4).
Resume
1. The distinctive feature of the presented approach to the Stefan problem modeling from the known methods consists in applying of the conjugate equation on the moving boundary instead of
traditionally used the boundary conditions of the IV type (the Stefan condition). In the new problem statement the heat exchange processes for the frozen and thawed zones and for the moving boundary are described by the same equations - the heat equations.
2. The problem numerical solution with the conjugate equation allows explicitly assign the frost penetration (thawing) moving boundary location with shock-capturing method implementation without the space-time grid reorganization on each new step on time.
3. The conjugate equation at each new time step gives the opportunity to define not only the moving boundary location, but the water aggregation state on the boundary.
4. Stage width of the frost penetration step graph gives pictorial presentation of occurring process dynamics, and for a concrete value gives the value of freezing time of a wet ground layer on the moving boundary.
5. Temperature-humidity compensation, brought in the problem algorithm at the point-to-point motion of the boundary, allows to comply the energy conservation law at the numerical solution level.
6. The numerical experiment results prove the Stefan problem solution method with the conjugate equation appliance doesn't yield in the calculation accuracy to known the Stefan problem solution methods in its classical formulation.
1 СНиП 2.02.01-83.* Основания зданий и сооружений.
2 СНиП 23-01-09. Строительная климатология.
3 СНиП 2.02.04-88. Основания и фундаменты на вечномерзлых грунтах.
References
[1] Лыков А.В. Тепломассообмен: Справочник. М.: Энергия, 1971. 560 с.
[2] Никитенко Н.И. Исследование нестационарных процессов тепло- и массообмена методом сеток. Киев: Наук. думка, 1971. 266 с.
[3] Никитенко Н.И. Сопряженные и обратные задачи тепломассопереноса. Киев: Наук. думка, 1988. 240 с.
[4] Основы мерзлотного прогноза при инженерно-геологических исследованиях / ред. В.А. Кудрявцев. М.: Изд-во МГУ, 1974. 431 с.
[5] Рыжиков А.А. Теоретические основы литейного производства. М.: Машгиз, 1954. 236 с.
[6] Самойлович Ю.А. // Металлургическая теплотехника, Свердловск: Среднеуральское кн. изд-во, 1965. Вып. 12. С. 114-137.
[7] Колесников А.Г. // ДАН СССР. 1952. 82. № 6. С. 889-891.
[8] Красовицкий Б.А., Шадрина А.П. // Теплофизика и механика материалов, природных сред и инженерных сооружений при низких температурах. Ч. 1 Теплофизика и механика природных сред и материалов. Якутск: Ин-т физико-технических проблем Севера СО РАН. 1974. С. 99-108.
[9] Богословский П.А., Битюрин А.К., Хорьков В.И. и др. // Инженерное мерзлотоведение в гидротехнике: Материалы конференций и совещаний по гидротехнике. Л.: Энергоатомиздат, 1989. С. 97-100.
[10] Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1972. 735 с.
[11] Crank J. Free and Moving Boundary Problems. Oxford: Clarendon Press. 1984. 425 p.
[12] Попов Ф.С. Вычислительные методы инженерной геокриологии. Новосибирск: Наука. Сибирская издательская фирма РАН, 1995. 136 с.
[13] Douglas J, Gallie G.M. II Duke Math. J. 1955. Vol. 22. № 4. P. 557-572.
[14] Васильев Ф.П. II ДАН СССР. 1964. Т. 157. № 6. С. 1280-1283.
[15] БудакБ.М., Васильев Ф.П., ЕгороваА.Т. II Вычислительные методы и программирование. М.: Изд-во Моск. ун-та, 1967. Вып. 6 С. 321-341.
[16] Javierre-Perez E. Literature Study: Numerical methods for solving Stefan problems. Delft: Delft University of Technology, Report 03-16, 2003. 94 p.
[17] КрасношлыкН.А., Богатырев А.О. II Вюник Черкаського ушверситету. Сер. Прикладна математика. 1нформатика, 2011. Вип. 194. С. 11-31.
[18] Будак Б.М., Гольдман Н.Л., Успенский А.Б. II Решения задач типа Стефана. М.: МГУ, 1972. Вып. 2. С. 3-23.
[19] Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.
[20] Самарский А.А, Моисеенко Б.Д. II Журн. вычисл. математики и мат. физики. 1965. Т. 5, № 5. С. 816-827.
[21] Будак Б.М., Соловьева Е.Н., Успенский А.Б. II Журн. вычисл. математики и мат. физики. 1965. Т. 5, № 5. С. 828-840.
[22] Вабищевич П.Н. Численные методы решения задач со свободной границей. М.: Изд-во Моск. ун-та, 1987. 164 с.
[23] ПоповЮ.А., Ращупкин Д.В., Пеняскин Т.И. Гидромеханизация в Северной строительно-климатической зоне. Л.: Стройиздат, 1982. 224 с.
[24] Попов Ю. А., Завалишина Т.В., Турантаев Г.Г. и др. II Изв. вузов. Строительство. 2004. № 10. С. 107-112.
[25] Завалишина Т.В. Энергосберегающая технология зимнего бетонирования строительных конструкций. Новосибирск: НГСУ (Сибстрин). 2003. 132 с.
[26] Мухетдинов Н.А. II Изв. ВНИИГ. 1971. Т. 96. С. 205-218.
[27] Богословский П.А., Ставровский А.П. II Гидротехническое строительство в районах Крайнего Севера: труды координационных совещаний по гидротехнике. Красноярск, 1975. Вып. 101. С. 90-93.
[28] Гольдин А.Л., Рассказов Л.А.Проектирование грунтовых плотин: Учеб. пособие М.: Энергоатомиздат, 1987. 311 с.
[29] Клейн И.С. II Инженерное мерзлотоведение в гидротехнике: Материалы конференций и совещаний по гидротехнике. Л.: Энергоатомиздат. Ленинградское отд-ние, 1989. С. 97-100.
[30] Арсеньева А.П., Богословский П.А., Февралёв А.В., Янченко А.В. II Инженерное мерзлотоведение в гидротехнике: Материалы конференций и совещаний по гидротехнике. Л.: Энергоатомиздат. Ленинградское отд-ние, 1989. С. 109-111.
[31] Битюрин А.К., Соболь С.В. II Гидротехническое строительство. 1993. № 11. С. 22-25.
[32] Битюрин А.К. II Известия вузов. Строительство. 1997. № 8. С. 53-58.
[33] Хаглеев Е.П., Хаглеев Ф.П. II Гидродинамика больших скоростей. Красноярск: КрПИ, 1992. С. 63-69.
[34] Хаглеев Е.П. // Вестник КГТУ Гидродинамикабольшихскоростей (теплоэнергетика). Красноярск: КГТУ, 1996. Вып. 3. С. 63-75.
[35] Хаглеев Е.П., Хаглеев П.Е. // Вычислительная математика, дифференциальные уравнения, информационные технологии: материалы Междунар. конф. Улан-Удэ, 2009. С. 244-252.
