CC BY
38
Вычислительные технологии
Том 9, № 5, 2004
NUMERICAL MODELLING OF CONVECTION OF ISOTHERMALLY INCOMPRESSIBLE FLUID UNDER LOW GRAVITY IN DOMAIN WITH FREE
BOUNDARY*
O.N. GONCHAROVA Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk and Altai State University, Barnaul, Russia e-mail: [email protected]
Численно исследуется тепловая гравитационная конвекция в жидкостях в условиях микрогравитации. В качестве математических моделей конвективных движений используются классическая модель Обербека — Буссинеска и альтернативная модель микроконвекции, характеризующаяся свойством несоленоидальности поля скоростей. Стационарная гравитационно-термокапиллярная конвекция рассматривается в полукруге со свободной границей. Проводится сравнение численных результатов исследования. Выявляются качественные различия в топологии течений в случае, когда граничный тепловой режим имеет локальную особенность. Численное исследование проводится для различных значений чисел Прандтля, Марангони и Рэлея.
Introduction
We consider the mathematical models of the thermal gravitational convection of liquids assuming smallness of the microconvection parameter. They are the classical Oberbeck — Boussinesq model of convection [1] and the model of microconvection of isothermally incompressible liquid known since 1991 (see [2, 3]). It was noted in [2, 3] that the approximation of Oberbeck — Boussinesq is not valid for a description of convection, if the parameter of microconvection is rather small. The microconvection parameter characterizes a ratio of the velocity orders produced by liquid expansion and buoyancy factor. It is equal to n = gl3/vx, where l is a characteristic linear scale of a region occupied by liquid, v, x are the coefficients of kinematic viscosity and thermal diffusivity, g = \g\, g is a gravity acceleration. The term "microconvection" was introduced in [2] to characterize the fluid flows with small values of the parameter n and therefore to describe the convective fluid flows under low gravity, in microscales or in liquids with large product of the coefficients of viscosity and thermal diffusivity. The microconvection parameter n is an additional criterium of similarity regarding to the Rayleigh Ra and Prandtl Pr numbers. The next specific term as "isothermally incompressible liquid" was used to determine
*The research was carried out in the frame of the grant 902.2003.01 of the President of the Russian Federation for support of the leading scientific schools.
( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
the liquid with the state equation, when its density p depends on temperature T only (see [2-4]).
Deriving both models of convection we proceed from the exact laws of conservation of mass, impulse and energy [1, 3]. The Oberbeck — Boussinesq system of the equations is a result of simplification of the complete conservation laws due to hypotheses: the density of liquid depends linearly on its temperature p = po(1 — ftT); motion is similar to motion of incompressible liquid and the velocity field V is considered to be solenoidal; in the equation of impulse the density variation is taken into consideration approximately. Additionally, the contribution of dissipation function and pressure forces are considered to be negligible and all transfer coefficients are constant. Here ft is a thermal expansion coefficient.
Then the Oberbeck — Boussinesq equations of the thermal gravitational convection can be written in non-dimensional form as follows:
divV = 0; (1)
^ ^ ^ 1 Ra
Ft + V ■W = —Vp' + - A V — R^pr 9oT ; (2)
Tt + V ■ VT = —^-AT. (3)
RePr
n
Here p is used for pressure, p' is modified pressure, such that p' = p — r 2p g0 ■ x, Re is the
Reynolds number, g0 = g/g.
The model of microconvection is based on the exact laws of conservation of mass and impulse [3]. The equation of energy is simplified due to hypothesis about neglecting of the contribution of dissipation function and pressure forces. All transfer coefficients are considered again to be constant. The alternative model is characterized now by a non-solenoidal velocity field V. Using dependence of liquid density on its temperature of type p = p0/(1 + ftT) the system of equations can be rewritten in form, when the modified velocity vector W becomes solenoidal [2, 3]. We introduce the system of the microconvection equations in non-dimensional form as follows:
divW = 0; (4)
-*-*-* F -* -* R2
Wt + W ■ VW + -^z(VT ■ VW — W ■ VT) + ^ 2 2 (AT ■ VT — V|VT|2) =
= (1+eT K-v?+f1^AWV')- R?k 9«T ; (5)
£ 1 Tt + V ■ VT + R^lVT|2 = (1 + £T) R^AT- (6)
v' 1 \ £
Here q is modified pressure ( q = p' — (1---p^j r AT, v' is a second viscosity) and
W = V—ReprVT •
The property of solenoidality of the modified velocity makes it possible to introduce an analog of stream function for two-dimensional and axis-symmetrical problems and to carry out the calculations of the convective flows in the variables "(modified) stream function — (physical) vorticity".
From physical point of view both above mentioned dependences of density on temperature are practically equivalent. In real convective flows the maximum values of ft|T| do not exceed
10-2 [1, 3]. The interest to the alternative models of convection has grown recently. These models are to play their role, for example, in the detection of so-called non-Boussinesq effects. The explanation of some experiments made on the orbital stations find no confirmation in the calculations using classical mathematical modelling.
At first the model of microconvection was proposed and used for the analytical and numerical research of non-stationary convection in the closed domains under low gravity [2, 5, 6]. Quantitative and qualitative differences in the characteristics of the non-stationary flows computed with classical and alternative models are confirmed. Taking into account the non-solenoidality for the stationary problems in the closed domains leads to the corrections of the order of Boussinesq number [7]. We continue the investigations of the microconvection problems in the domains with free boundaries in case, when the boundary thermal regimes have local singularity, and model additionally fast changing boundary regime (see [8, 9]). The calculations are carried out for the different Prandtl, Marangoni and Rayleigh numbers. The topology of convective flows is investigated and those situations, when the topology of flow can differ, are presented.
We will consider the stationary gravitational-thermocapillary convection in a semicircular domain with free flat boundary. Under conditions of low gravity and in case, when the parameter responsible for deformation of free surface by thermo-capillary forces (the capillary number) is rather small, a non-deformed free boundary is considered. It can be approximately defined as a boundary of capillary balance. The correction to the free boundary can be found from dynamic condition on free boundary.
Remark 1. In this paper the following usual notations are used: r for a radial coordinate, for an angular coordinate, и for a vorticity, ф for a stream function or for modified stream
1 дф дф function. Then v = -7т— is the radial component of velocity, u = —-p— is the circumferential
r d^ dr
component of velocity. We introduce the definitions of the non-dimensional parameters: Ra =
en, ^ = в7*, Re = v*//v, Pr = v/x, where e is called the Boussinesq number and "star" is
used to notice the characteristic values of functions. A choice of the characteristic values can
be done according to [10].
1. Formulation of problems
Stationary gravitational-thermocapillary convection is investigated in the semicircle of type
0 < r < R < п < < 2п.
The diameter of semicircle (^ = п, ^ = 2п, 0 < r < R) is free boundary and the semi-circumference (r = R, п < < 2п) is rigid boundary with given heat flux through this boundary. The equations of convection for both mathematical models (1)-(3) and (4)-(6) can be rewritten in the variables ф — и.
1.1. Classical Oberbeck — Boussinesq model in terms ф — и
The equations (1)-(3) considered in stationary case are written in the polar coordinates:
Ди Re (vdu иди^ Ra (dT cos IdT sin ^ =0 (7)
U \V dr r д<£ / Ma \ dr COS ^ rd^ S j ,
Дф + и = 0, (8)
a^ ™ « dT udT\ AT — Ma( v—- + 1 =0, dr r
(9)
where Ma = RePr is the Marangoni number.
To realize the fast changing temperature regimes the local singularity of thermal flux through free boundary is created. In this connection two types of boundary conditions for both models are studied. Let us indicate them symbolically: Variant I or the basic variant with no "splashes" on free boundary and Variant II or the variant with "splash" on free boundary. "Splash" means an action of local singularity for the boundary regimes. Additionally the different changes of the temperature regimes on rigid boundary are considered.
We introduce the boundary conditions for temperature on the rigid boundary as follows:
dT
— = Tg cos 7<£, y = {1, 2, 4}.
Variant I. The boundary conditions on the free boundary 0 < r < R, ^ = n, ^ are considered for stream function and vorticity
^ = 0, u
ff, P = 2n, dr ' r '
dT =
(10) = 2n
(11)
and for temperature
£ = 0.
On the rigid boundary r = R we obtain from the no-slip conditions as usually that
*=o, dr=0.
(12)
(13)
Variant II. For the second variant we consider (11) and write the boundary conditions for temperature on the free boundary 0 < r < R, ^ = n, ^ = 2n as follows:
ÔT
0, = n, ^ = 2n (r = R*),
R*Tb, = 2n (r = R*).
(14)
On the rigid boundary r = R we keep the conditions (10), (13).
Remark 2. For the Variant I we have the problem statement consisting of the equations (7)-(9) and boundary conditions (10)-(13). For the Variant II we consider the equations (7)-^) and the boundary conditions (10), (11), (13), (14).
1.2. Model of microconvection in terms ^ — w
The non-dimensional equations (4)-(6) for stationary problem in the polar coordinates can be written as follows:
rrr, A ^ ( dw wdw\
[1 + fT] Aw — Re v— + +
\ dr r d^y
1 dTd? I
r dr r dr
dT / u\ 1 dT /A v — Au — — — -—- Av — —
dr
r
+
r
r
Ra /dT 1 dT -7— cos --7— sin
r d^
Ma V dr
fuAT —du 1 dldu^
Pr y dr dr r2 d^ d^ y
e2
MaPr
1dT dAT 1dT dAT -----1----
r dr d^ r d^ dr
0;
A" + u = 0;
[1 + eT] AT - Ma
dT «dT
dr r d^
- el V TI
0.
(15)
(16) (17)
Here q = Re q. In the equations (7)-(9) and (15)-(17) we use following notations for the differential operators in the polar coordinates
A 1 d d 1 d2 A =--r--1---,
r dr dr r2 d^2
V =| -d, I—
dr r d^
We introduce again the boundary conditions (10) on the rigid boundary.
Variant I. For the first variant we obtain the similar to (11), (12) boundary conditions on the free boundary 0 < r < R, ^ = n, ^ = 2n for the modified stream function, vorticity and for temperature. On the rigid boundary r = R we have now for the modified stream function that
e m 1 . d^ 1 £ dT
" = — R^i-Tq— sin ——
Ma
Y
dr R Ma d^
(18)
Variant II. For the second variant on the free boundary 0 < r < R, ^ = n, ^ = 2n we consider the conditions (11) for the modified stream function and vorticity and the condition
(14) for temperature. On the rigid boundary r = R we consider the conditions of type (10), (18). Remark 3. The problem statement for the Variant I consists of the equations (15)-(17)
and boundary conditions (10), (11), (12), (18). For the Variant II we consider the equations
(15)-(17) and the boundary conditions (10), (11), (14), (18).
4
2
2. Numerical procedure
The mathematical models of convection described in Section 2 are investigated numerically. Numerical procedure for their solution is carried out with help of the finite difference scheme based on the alternating direction method. This method formally has second approximation order [11]. The convective terms taken from the proceeding iteration layers are approximated against flow. Actually, the scheme of first order is obtained. The method offered for investigations was approved on the test problems and calculations of non-stationary convection in the circular domains [12] (see also [8, 9]).
For the equations (7), (9) or (15), (17) the scheme of computation is written in the following form:
U fc+1/2 — TT k
U 05- U = au [AxUk + A2Ufc+1/2] + auFk,
Uk+1 _ Uk+1/2
u U - = au [A1Uk+x +A2Uk+1/2] + auFk, (19)
where U = , Uk = U (tk), A1 and A2 are the difference operators, which approximate the differential operators in A , Au = Au for the Oberbeck — Boussinesq model, Au = Au(1 + eTk)
0.5t
for the model of microconvection. Here A^ in an iteration parameter, and Fk includes all terms in the left-hand sides of the equations (7), (9), (15), (17) beginning from the second one and calculated on the previous layer.
Numerical solution of the equations (8) or (16) at each step = kr, k = 1, 2,..., is obtained by similar iteration scheme
- " = A^ (A1^s+1/2 + A2/ + ,
0.5t
A^ (Ai^s+1/2 + A2^s+1 + w^1) (20)
0.5t
with another iteration parameter A^.
To realize the stated above calculation scheme a difference grid is introduced:
r„ = (n - 1)h (n = 1,..., N +1), h = R/N,
= (m — 1)a (m = m,..., M + 1), a = 2n/M (ma = n). Here we use the notations f (rn, = fn,m and A1f , A2f are the difference operators:
Aif = -1
rn
fn+1,m fn,m fn,m fn- 1,m
' n+1/2--h2 'n-1/2 ^2
+ fn
A 2 f
2 J =
r2 a2
n
The idea of approximation against flow is used for the approximation of the convective terms:
d J u d J w— +----
dr r d^
W fn+1,m fn-1,m + un,m fn,m+1 fn,m-1
n'm 2h rn 2a
+
+ fn + fn
+1 vn,m1 rT I o •
2h rn 2a
The first derivatives on boundary are approximated by one-side differences.
In order to determine the boundary conditions for vortex on the rigid boundary the conditions of Thom type [13, 14] are introduced with help of the Taylor expansions and Poisson equations (8), (16) considered on the boundary:
^N +1,m = — h22 "n,™
2 / £ dT f 1 + M e T • (1 +2R 1\
WN= — h^— Ma d^R + hRj— Ma Tg sin 7M RY+^ y) •
These conditions are written for classical model and for the model of microconvection respectively.
We present a general scheme for solution of the problems consisting in the realization of the following stages:
1. External iteration process consists in successive calculations of the functions Tfc+1, wfc+1 from the equations (7), (9), (15), (17). Moreover the Thomas algorithm in the direction ^ is realized on every intermediate (k + 1/2) layer. On the basic (k + 1) layer the Thomas
r^j
algorithm is realized in the direction r. The initial data are determined by the rest state of type T := T0 = const, u := 0, ^ := 0.
2. Internal iteration process of calculation of from the equations (8), (16) is introduced on every (k + 1) iteration layer with the alternated sequence of the Thomas algorithms. After an end of the iterations s = S it is considered that with a given precision e^ the values of ^ are determined on the (k + 1) layer, such that ^fc+1 = .
The iteration processes are considered to be convergent, if the criteria of convergence are fulfilled:
max ¡/n+m - J < ef max ¡Z^
n,m ' ' n,m '
where i is an iteration number, ef is a precision of the calculations of /i+1 (see [13, 14]). We use an additional examination for the fulfillment of the boundary conditions [14] with help of
e = max |un) - un+i,m(^Nrm)|
m '
and a condition that the stationary flow is considered to be achieved, if no less than K of the external iterations are fulfilled. An achievement of stationary solution is rather delicate question. Therefore an examination of an exit on the stationary regime with help of a perturbation of "an initial" approximation was used. Returning to the earlier achieved "stationary" state was observed for all presented cases.
Remark 4. The questions connected with a correction of the free boundary and with a stability of the flows can be solved according to [10, 15]. If H(x) is a deviation of free boundary from the position y = 0, — R < x < R, an equation for this correction can be written as follows:
R
iP—iie=— caRH" H(±R) = °. /Hd* = °-
-R
Here "prime" denotes a derivative on x, v2 is expressed by radial and tangential velocity components v, u, is a deviation of the pressure from the balanced level, a is a surface tension and ca is the capillary number.
It should be noted that the stability of convection flows in domains with free boundaries is of great interest (e.g. [15]). Computational results presented here show that the instability phenomena will occur at higher values of Re and Ma that were reached in this paper.
3. Numerical results
The calculations are performed for the physical liquids similar to glycerin and melts of silicon and called symbolically Glycl, G/yc3, Si/. The main parameters of the substances can be found in the Table. The solutions are computed on the grids with 41 x 41, 81 x 81 and 161 x 161 mesh points. The results of calculations are shown in Fig. 1-3. The value of radius is R =1 (cm). The values of the parameters in the boundary conditions are: TB = 70 and TB = 150, TG = 35, To = 35, R* = 0.45.
We describe now the results of numerics.
Variant 1 (Basic variant). Case 7 = 1. In this case the computations for two alternative models show only some quantitative differences and practically the same qualitative pictures. The stationary solution has the one-vortex structure for the liquids G/yc3, Si/ and the two-vortex flow of type "two small vortices in one" for the liquid G/yc1.
Parameters of the problems
Pr Ma Re Ra n e
Glycl 104 3 ■ 102 3■10-2 1.5 ■ 10-3 10-1 1.2 ■ 10-2
Glyc3 104 1 10-4 1.5 ■ 10-3 10-1 1.5 ■ 10-2
Sil 4■10-3 1 2.5 ■ 102 2■10-4 1 2■10-4
In the basic case with 7 =1 a stability and an experimental order of convergence r of solution of the difference problem was tested due to the Runge rule (see [16, 17]). To estimate the error of the numerical results we can consider some measurement ri, r2, r3 on subsequently refined grids, here (i = 1) : 41 x 41, (i = 2) : 81 x 81, (i = 3) : 161 x 161. We determine the quantities r», which characterize so-called motion intensity max These motion characteristics are calculated
for the liquid G/yc1: r1 = 0.0345, r2 = 0.0303, r3 = 0.0291. We obtain the experimental order of
convergence r ~ 1.8 calculated as follows: r = ln(|r2 — r1 |/|r3 — r21)/ln2. The estimated relative error of the motion intensity is about 5%. (For its calculation we have used [1/(1 —
(1/2)r)] ■ [|rs — r2|/rs].)
In the Case with 7 = 2 the calculations made for the liquids G/yc1, G/yc3, Si/ show the two-vortex structures of flow. The same pictures of topology of flow and of temperature field are observed in results due to both mathematical models.
Case 7 = 4. In figures (Fig. 1, a — c) the topology of flows and the temperature fields are presented. By the calculations made for the liquid of type G/yc1 the four-vortex structure of flow is observed (see Fig. 1, a). In Fig. 1, c the two-vortex structure of flow for the liquid of type Si/ is presented. The typical group of the isotherms is given in Fig. 1, b. The orders of non-dimensional velocities are ~
10-2 — 10-1 inside of region and ~ 10 1 at the free boundary
c
Fig. 1. Variant 1, 7 = 4: a — topology of flow for Glycl; b — temperature field for Glyc1; c — topology of flow for Sil.
for G/yc1, ~ 10-1 — 10° inside of region and ~ 10° at the free boundary for Si/.
So, for basic Variant I the qualitative differences in topology of flows computed by different mathematical models are not observed.
Variant 2. (Additional modelling of fast changing temperature field by creating local singularity of the thermal flux on free boundary.)
Case y =1. We consider rather weak singularity of Gauss type with TB = 70. The differences in the computations appear in the investigations of flows for the different liquids. For Si/ we obtain the two-vortex flow structure of type "two vortices in one", where an inner small vortex is located by the singularity point (see Fig. 2, a). For G/yc3 both mathematical models give the one-vortex structure of flow (Fig. 2, b) The temperature field computed for G/yc3 is presented in Fig. 2, c. The similar temperature picture is also obtained for Si/. The temperature variation of G/yc3 is in range from 21 to 50 for the Oberbeck — Boussinesq model and in range from 18.5 to 55 for the microconvection model. The temperature variation of Si/ is in the interval [25; 45] for the Oberbeck — Boussinesq model and in the interval [21; 50] for the microconvection model. The orders of the non-dimensional velocities are ~ 10-1 — 101 inside of region and ~ 101 at the free boundary for G/yc3. The orders of velocities are ~ 10-1 — 10° inside of region and ~ 10° at the free boundary for Si/.
Case y = 4. The more distinct differences in the calculations using the alternative models are appeared for the liquid of type Si/ for the singularity with TB = 150. In Fig. 3, a the topology of flow calculated according to the model of microconvection is shown. The rather complicated two-vortex flow structure can be shown. The right big vortex is the vortex of type "two vortices in one". In Fig. 3, c the topology of flow computed by the Oberbeck — Boussinesq model is shown. In this case the flow is characterized by two vortices. They are of different sizes and have approximately equal intensity. The isotherms field is presented in Fig. 3, b. One can speak about similar pictures for the temperature fields. The temperature variation is in range from 26 to 45 for the Oberbeck — Boussinesq model and in range from 26 to 49 for the
Fig. 2. Variant 2, y =1: a — topology of flow for Sil; b — topology of flow for Glyc3; c — temperature field for Glyc3.
a b
c
Fig. 3. Variant 2, 7 = 4: a — model of microconvection (topology of flow for Sil); b — temperature field for Sil; c — Oberbeck — Boussinesq model (topology of flow for Sil).
microconvection model. The orders of non-dimensional velocities are ~ 10 2 — 100 inside of region and ~ 100 at the free boundary.
The quantitative results obtained by both models are quite close to each other. The values of velocities differ approximately on 15%.
Conclusions
Our mathematical modelling of the convection can be called alternative approach in the convection theory. The most bright qualitative and quantitative differences from the classical results or the non-Boussinesq flow effects are obtained by the simulations for the non-stationary problems of microconvection [2, 3, 5, 6]. According to analytical results of V.V. Pukhnachov [7] the non-solenoidality for the stationary problems of microconvection in the closed domains leads to the corrections of the orders of Boussinesq number.
The purpose of the paper is not only to demonstrate the topology of flows for different liquids in domain with free boundary. We want to show numerically the possibilities for creating differences from the classical theory in stationary case. The qualitative differences in flow characteristics for stationary problems with free boundary in a semicircular domain can be observed, when boundary thermal regimes have local singularity. Modelling of local singularity should be combined with fast changing heat flux through fixed boundary (for instance, by condition of type (10) with 7 > 1). The differences in the quantitative flow characteristics obtained by different mathematical models will not be so essential as in the non-stationary problems.
References
[1] Joseph D. Stability of Fluid Motions. Berlin, Heidelberg, N.Y.: Springer-Verlag, 1976. 638 p.
[2] PüKHNACHOV V.V. Microkonvection in vertical layer // Proc. of Russ. Acad. of Sci. Fluid and Gas Mechanics. 1994. N 5. P. 76-84.
[3] Applications of Group-theoretical Methods in Hydrodinamics / V.K. Andreev, O.V. Kaptsov, V.V. Pukhnachov, A.A. Rodionov. Dordrecht, Boston, L.: Kluwer Acad. Publ., 1998. 408 p.
[4] Nadolin K.A. Convection in horizontal layer under the inversion of specific volume // Proc. of Russ. Acad. of Sci. Fluid and Gas Mechanics. 1989. N 1. P. 43-49.
[5] GONCHAROVA O.N. Microconvection in weak force fields. A numerical comparison of two models // J. of Appl. Mech. and Technical Phys. 1997. Vol. 38, N 2. P. 219-223.
[6] GONCHAROVA O.N. Simulation of microconvection in a long rectangle // J. of Computational Technologies. 2000. Vol. 5, N 5. P. 26-37.
[7] PuKHNACHOV V.V. Stationary problem of microconvection // J. of Continuous Medium Dynamics. 1996. N 111. P. 109-116.
[8] GONCHAROVA O.N. Numerical simulation of microconvection in domains with free boundaries // J. of Appl. Mech. and Technical Phys. 1997. Vol. 38, N 3. P. 386-390.
[9] GONCHAROVA O.N. Microconvection in domain with free boundary // J. of Computational Technologies. 2000. Vol. 5, N 2. P. 14-25.
[10] Zebib A., HomskyG.M., Meiburg E. High Marangoni number convection in a square cavity // Phys. Fluids. 1985. Vol. 28, N 12. P. 3467-3476.
[11] Samarsky A.A., Vabischevich P.N., Matus P.P. Differences Schemes With Operator Multipliers. Minsk: 1998. 442 p.
[12] Voevodin A.F., Goncharova O.N. Calculation of free convection with variable gravity field // J. of Continuous Medium Dynamics. 1984. Vol. 67. P. 21-28.
[13] Roache P.J. Computational Fluid Dynamics. Albuquerque: Hermosa Publ, 1976. 616 p.
[14] Tarunin E.L. Computational Experiment in the Problems of Free Convection. Irkutsk: Irkutsk State Univ., 1990. 228 p.
[15] Myshkis A.D., Babsky V.G., Kopachevsky N.D. et al. Low-gravity Fluid Mechanics. Mathematical Theory of Capillary Phenomena. Berlin: Springer-Verlag, 1987. 504 p.
[16] Ostapenko V.V. Finite diffenerce scheme of high order of convergence at a nonstationary shock wave // Siberian J. of Numer. Mathematics. 1999. Vol. 2. N 1. P. 47-56.
[17] Dörfler W., Goncharova O., Kroner D. Fluid flow with dynamic contact angle: Numerical simulation // ZAMM. 2002. Vol. 82. N 3. P. 167-176.
Received for publication February 16, 2004
