- Email: [email protected]
A numerical simulation scheme for subsonic flows in a viscous compressible gas
U.S.S.R. Comput.Maths.Math.Phys.Vol.28,No.6,pp.l77-182,lg88 Printed in Great Britain
0041-5553/88 $~O.OG+O.OO 01990 Pergamon Press plc
A NUMEhICAL SIf:lULATION SCHEKE FOfiSUBSONIC FLOWS IN A VISCOUSCOMPRESSIBLE GAS*
V.A. GONCHAROV,
V.N. KRIVTSOV
and A.A. CHAFAKHCH'YAN
A difference scheme is proposed for two-dimensional flows of a compressible gas, comparable in efficiency with schemes for the equations of an incompressible fluid. 'Ihe efficiency of the scheme is determined by an original iterative process, using a procedure for eliminating for the equations of an incomthe pressure, which is characteristic A comparison with other schemes is made for a model pressible fluid. convection problem. When one is studying subsonic time-dependent flows of a gas under conditions of large temperature drops, the incompressible fluid approximation may prove unsatisfactory. Such flows are described by the full Navier-Stokes equations. It is natural to require that a difference scheme for these equations perform well even in the limiting case of practically incompressible flow - at least as well as the standard difference schemes for the equations of an incompressible fluid. This means that: 1) the volume of computational work per time step should be comparable, and 2) computations at the same time step should achieve the same The scheme proposed below meets these requirements. accuracy. In difference schemes for the incompressible fluid equations, each time step as a rule involves the solution of a difference analogue of Poisson's equation. Among the schemes available for the equations of a compressible gas are those that completely split the process with respect to the space variables (see, e.g., /l/j. Such schemes are even more efficient per time step than schemes for the equations of an incompressible fluid. However, it turns This is due to the need to calculate with a out that the second condition is not satisfied. time step-size T-h/u, where h is the grid step-size and u is the velocity of the gas. Corresponding to this step-size, however, one has large values of the Courant number (rclh, where c is the velocity of sound), which implies an increase in that part of the approximation This was shown graphically in error introduced into the scheme by the splitting procedure. /l/, where computation with a time step-size 16 times greater than that required by the condition produced a very distorted solution. Courant The difficulties involved in devising an efficient scheme for the subsonic flows of a compressible gas within the limits of the full Navier-Stokes equations have inspired the development of an approximate theory (see e.g., /2, 3/j, which permits a simpler construction This approach has undoubtedly been of difference schemes with the required properties. successful when applied to various problems. Nevertheless, we do not believe that this precludes the development of new computation schemes for the full equations for a compressible flows of a gas. Suffice it to mention that the application of the theory of time-dependent heat-conducting gas has not been adequately substantiated, as the pressure variation in such flows need not be of the same order of magnitude as the square of the Mach number.
1. We shall consider the havier-Stokes equations for a heat-conducting gas, including the effect of the force of gravity with a constant acceleration g. The system is completed by the equation of state of an ideal gas P=RpT, where P is the pressure, p is the density, The equations will be presented in a form T is the temperature and R is the gas constant. suitable for the convection problem. Let Z,pO, T0 be the characteristic length of the region, Together with the total pressure P we introduce the density and temperature, respectively. where Pst is the hydrostatic pressure of a gas of density excess pressure p=P--RpoTo-Pst, flow velocity is defined by the formula V=(gZ)"'. We introduce nonPO. The characteristic dimensional variables r&r’,
u=l.u’,
t=t’Z/V,
T=T,T’,
(P, p, P&=
poV2(P’, P',r&L Here r is the radius-vector of a point of the space, u is the velocity vector and t is the variables are indicated by primes, which will be omitted henceforth. time; non-dimensional The Reynolds, Prandtl and Mach numbers are defined by the formulae Re=p,I'Z/y,
Pr=l*C,lh,
MZ= V*/~I?T,,
y=c,/c”,
where ).Iand A are the viscosity and thermal conductivity, respectively, which will be assumed constant for simplicity,and C, and C. are the specific heats at constant pressure and constant volume. *zh.vychisl.~at.mat.Fiz.,28,12.1858-1866.1988 177
178
The transport operator in the equations can be reduced to various forms by identical For the energy equation we shall use the divergent form of the transport transformations. operator, and for the equations of motion, the skew-symmetric form (see /4/). This necessitates w=p"u. In the equations presented below the introduction of a new unknown vector-function we have used both variables u and W;this is dictated by the structure of the difference scheme 11-p'". After these transformations, We shall use the abbreviated notation proposed below. the equations are written as follows:
q $i+K(u,w) C
=-gradp+(p--l)l+ I
A(m/n)-i-$graddiv(w/n)
P-p+p,
$+div(r~w)=O, Do=-=
Irz
M*y(r-l),
C,To
1,
t+D-DpT,
D,=&,
D=+=(~M’)-I.
Here 1 is the unit vector in the direction of the gravitational force, Q, is the work of the is a skew-symmetric operator, for which we have viscous forces, K(u, w)=(ugrad)w+0.5mdivu This property of the operator K will also be carried over the identity wK(u,w)=div(uwz/2). to its difference approixmation (see /4/j, so that the schemes obtained for the equations of motion will involve a transport operator that conserves kinetic energy. The boundary conditons are as follows: for the velocity - no-slip conditions, and for the temperature - a condition of the first kind on certain sections of the boundary and thermal insulation on the others. the 2. We shall consider two-dimensional plane or axially symmetric flows, with P Differences between these two cases will be radial coordinate and s the axial coordinate. indicated in the formulae by the index v:v=O for plane flow and 11=1 for axially symmetric In these variables the differential operators are flow. grad =
alar IIa/a2 II1
divw=1dr'w,+aw, TV ar az ’
AT+-rrv$+$)T,
To fix ideas we shall assume that the force The flow region is a rectangle in the r, z plane. In the I‘, s plane we introduce a rectangular grid of gravity acts along the s coordinate. with constant step-sizes h, and h,, whose points are numered by integer indices i, j. Following /5/, we define the thermodynamic variables at the centres of the grid cells:(P, T, P)<+~,i+~h~ and the components of the vector W at the midpoints of the appropriate sides: WV,,. i+'br zaz, ‘+‘h. J. The values of the difference functions at other points of the grid are obtained by Linear interpolation from the nearest grid-points at which these functions are defined. The boundary of the region consists of several grid lines, on which the normal velocity To specify the temperature component is defined, as well as the thermal flux, if necessary. and the tangential component of the velocity vector on the boundary, we introduce fictitious cells, whose data are defined in such a way that linear interpolation together with the data in neighbouring cells gives the prescribed boundary value. In constructing the difference scheme, we shall use the idea of splitting wherever we consider it justified. To obtain second-order time accuracy, we will use the symmetric splitting scheme of /6/. The sequence of computations per time step is determined by the following difference formulae: (w,),+0.5K,h(U,, wr"")=O, (pT),f0.5
(w,),+0.5K,"(u,, f1'*05)=o,
div,"[u(pT)'.']=O,
(w,),+O.5K,"(u,, w,o")=O,
Re-’
(PT),=D,A,“T’,~,
A,‘“( w,O “/a),
(p),=O;
(lb)
(w,),i-0.5K.h(&, w,"~)=O,
(pT),+0.5div,"[z~(pT)~~j=O, (tc,),=O.5q-’
(p),=O;
(P) r=O; (w,),=O.5q-’
(la)
(2a) (P)
Re-’ A,*( to,” “/q).
(3rd (3b)
179
,
(2~),=0.5q-~Re-'A," (pT),=D,A,"T",
(~,).=0.511-'Re-'A,L(~),
q
(p).=O;
(4b)
(w,),=-$~rad,"div'(~),
(,?a)
( 1 1;
grad,"div"E
(wz)T=&
(5b)
$=[p(pO.J,T'.")-l]Z,,
n(w),=-grad"p"+$,
(44
(pT)~=--D,~0.5diP(~a5/7])f~yh,
(p)z+div"(qw)=O,
(6a) (6b)
P=p+P&.&D=DpT. Computations using formulae (l)-(5) are then carried out in the reverse order. Here (v)z=(qn(P~)/T, where T is the time step-size, 'pn is the "new" value of the difference function, to be determined at the current step of the splitting scheme, and 'p. is the "old" value, known from the previous step; @."=0.5(cpnt~o). The function p(p, T) at step (6a) is determined from the equation of state. The difference functions u and n are computed from the data on the previous time layer and not changed on the current time layer. The superscript h denotes difference differentiation operators with respect to the appropriate variable, the latter being indicated by subscripts r or z. These operators have a natural form for the grid: Mlh(&f)
[&‘“fl
lIJ+‘b=
[ Vu) *+‘/,fi+*(Jw I-*k-1 lj+1,, I Zh,r,‘l
iW)r-(rvl)r-l ,,,+'I>= -h,= I I-‘ ' 3+ !- L, 1,+'I* 1
(Pf)S+l-(rVl)
[rdl (fn++-f,+~IJ) -riy (f,+~,,-f,_s,JI, [A:fl z+a/>,, = r:+,,xh.* (A.“f)i+v,,,=
f+%[f~+I-2fj+fj-i1 h.’
’
f[fj+a,:-2f,+a+fi-al (A.“f)t.,+‘,:= It,* ’
km&” /Ir+‘,>.j =
fi+‘b.j+lh-f,+‘l.,l-‘h~ hz ’
f f+B,~+l-ff+‘lr., (div,”f) f+‘hJ+% = 11, ’ div” w* = div,"LU,+ div,"w,. ( I(', 1 In some of these formulae subscripts common to several variables have been written outside the brackets (the first subscript before the brackets, the second after them). The expression for the work of the viscous forces also has a natural form, not shown here (in convection problems this term is small compared with the others). The approximation of the operators K(tr,U?) was constructed as in /4/. Steps (l)-(4) are implemented using one-dimensional pivotal condensations. Step (5) is accomplished by a version of the method of alternating directions. The function wI in Eq. (5a) is assumed to be known. Then ZL',is found by one-dimensional pivotal condensations and substituted into Eq.(Sb) , which is solved by one-dimensional pivotal condensations for mz, and so on. for In practice it has been found that two or three iterations are sufficient convergence. We will now consider the system of algebraic equations arising at step (6). The unknowns 10,p, T,p. We use a well-known device /7/, in which the here are the difference functions equations are divided into two groups: the first includes the energy equation, and the second, the remaining equations. First, assuming T to be known, one solves the equations of the
180 second group and finds the other unknown functions. Then T is determined from the energy equation, and so on. In performing these computations it was assumed that the maximumadmissrelative change in temperature in one extraneous iteration was ible 6=1o-s; to achieve this it wassufficient to perform one or two extraneous iterations. The efficiency of the scheme is largely determined by the iterative process for solving the second group of equations, which will be consideredinSect.3. 3. We introduce a new difference vector function V'T)IlF. Then the system of equations of step (6), given the difference function T, is written in the form 2(v),=-gradhpO,'++,
(uL=(v-_gWo)l',
(7)
(P)~=-div"v,
(3)
p=DpT--D---Pst
(9)
The function $ in the equation of motion will be considered to be known. We first eliminate the difference functionspand p. To this end, we substitute (9) and substitute the result into (7). This gives 2(~),=0.5D~ grad"T div” v+q, (P-Q-0.5
gradh(po--Ps-,-0,5D
(8) into
(1Oa) gradh(po7’),
(lob)
where the subscript "o", as before, denotes functions known from previous steps of the splitting scheme. Rather than solve this system directly, we shall eliminate the pLessure from the equations of motion, a normal procedure when solving the equations of an incompressible fluid. We introduce the difference operator rothu=grad,'hu,-grad,'hu,, where grad."'and grad,'* are the operators of difference differentiation with respect to the respective variables, whose values are defined at the integer-indexed points of the grid: (grad,.hf),,
=
fy-fw,
(grad,.“f)
,j
=
h+;-h--y
I
It is readily Applying
seen that rothgradh=O. the operator rot*
to Eq.(7),
we obtain
rot"u=f,=lT rothIp+roth u. (11) Assuming now that (P)~ in (8) is given, we obtain a difference analogue, of the well-known problem of determining a vector field, given its divergence and curl. If the normal component of the vector on the boundary of the region is given, theproblem is not necessarily solvable for any function defining the divergence: the integral of the function must equal the integral of the normal component of the vector along the boundary. In the case of the difference equations (8) and (11) this property of the differential problem is represented by the fact that the number of equations (for a given function (p)* exceeds the number of unknowns by one. If the function (P)~ is properly chosen, the rank of the system of equations (8), (11) equals the number of unknowns, since the sum of all the equations of (8) becomes an identity, inthis case Zdi@u=O, where the summation extends over all cell centres and it is taken into consideration that the normal component of the velocity vanishes on the boundary. Below, (P)~) will always be caiculated from Eq.(8) for some given function V satisfying the boundary conditions, and therefore the system of Eqs.(8), (11) is consistent. We will now construct a procedure for eliminating one of the components of vfrom systen(8), (11) making use of the analogy with the familiar formula rotrot=graddiv-A for the differential operators. We formally extend the definition of Eq.(ll) to the boundary lines of the grid r=conSt, say by setting f,=O. To this end it is necessary to introduce fictitious cells near the boundaries and suitably define the function u, in them. Since v,=O on the boundary, we arrive at the condition VI,%,j=JJ., ‘i:, ,, (12) where the subscript '/*corresponds to the fictitious cells and '/z corresponds to the nearest interior cells. We will now consider the operator div?"', whose values are defined at the same points as the component v,,with the exception of the boundary lines z=const: (div,'" f) 1+%.,=
ff+l,jl~:+L--fUTiV “:+;, hr
.
It is readily seen that div,"'gradr*"=gradZhdivlh. Apply the operator operator grad," to ~q.(8) and then subtract one from the other. Poisson equation for v,: A."u,=f=-gradZh(p),-div,'"f,, A"=grad,"div,"+div,'h grad,'h.
div,'" to Eq.(ll) and the This gives a "difference" (13)
181 and a condition of type (12) on the lines The boundary conditions are U.-O equation was solved using the direct method of the ordinary Fourier transform
z=const. This /6/.
NU
4
b
a t
: 3
2 : 31 k 0
10
20 Fig.1
NU
4 t
a
Fig.2 We can now formulate the proposed iterative process for solving Eqs.(i')-(9). Fix some Substitute the resulting function function (pf=, say identically zero, and solve Eq.(13). uI into (loa), and then solve this equation for U, using one-dimensional pivotal condensations. pivotal Kow substitute the resulting function V, into (lob), and again use one-dimensional condensations to find a new v,. The last step of the iteration is to substitute the computed vector U into Eq.(8) and calculate a new function (p)%. In all our computations, the iterative process converged fairly rapidly - in two or three iterations - even for comparatively large Mach numbers (0.8 to 1.0). Counting also the extraneous iterations due to the change in temperature, one time step required solving Poisson's equation at most six times. As shown by comparison, the computation of strongly compressed flows takes all in all 2 to 2.5 times m3re time than the analogous computation with the Boussinesq approximation. 4. Th e proposed numerical method (method III) was compared with the explicit scheme of /6/ (method I) and the implicit splitting scheme of /l/ (method II) for the model problem of /l/. Consider the motion of a compressible gas in a closed square of unit length under the action of the force of gravity. The upper and lower walls are thermally insulated, but a constant temperature is maintained at r=l , T= 1.5 at on the side walls: T=i r=O. The velocity satisfies the no-slip condition at the boundary of the domain: u,=u,=o. The initial conditions are as follows: a linear temperature distribution in the gas at rest: T=l.S-0.5, uI,u,=O, p=D=const, corresponding to the steady state in the absence of gravity. bon-dimensional parameters: Re=250. Pr=O.71, y=1.4. The problem is solved for two values of the Mach number: M,z=1/320, MZ2=1/20. The result of the computation is the behaviour with time of the mean Nusselt number at the hot wall:
The computations were carred out with different time step-sizes on a 21x21 grid. To check the accuracy of the computations in the steady section, it was verified that Nuj,,o= NIlI,-,. In Fig.la the Kusselt number Nu(t) at the hot wall is plotted against time for M=M,, 1.=0.0056 and method II with r=2r., 42.,8% using method I with maximum admissible step-size (curve 2) and ~=32~. (curve 3). An eightfold increase in step-size (curve 11, r=16r. yields a deviation of less than 1.5% from the solution of the explicit scheme, and therefore A further increase in all these results are represented by the same curve on the graph.
182 step-size caused a marked change in the Nu(t) values (curves 2 and 3). Fig.lb shows similar olots of the Nusselt number Nu(t) obtained by method I with r=r.(curve (curve 2) , r=32r. (curve 3), z=64r. (curve 4) 1) and by method III with r=&., 4t., 8r., 16~. and r=128~. The solutions represented by curves 1, 2 and 3 differ by less than (curve 5). 2%. In computations with .M=M, using method III, the discrepancy between the mean Nusselt numbers of the cold and hot walls amounts to ANu/Nu=2%. This discrepancy determines the accuracy of the computations. Thus, the solutions corresponding to curves 1, 2 and 3 are identical, to within the accuracy of the computations. The discrepancy between the solution to -5%. The larger with ~=64r. and the solution obtained by the explicit scheme amunts difference between these solutions is obvious1 y due to the fact that the steady solutions of the splitting scheme depend on the time step-size 'c. It is also worth noting that all the curves obtained by method III duplicate the shape of curve 1, whereas the curves produced by method II show considerable fluctuations about the solution obtained from the explicit scheme. Similar results were obtained for the other value of the Mach number, M=Mt. Fig.2a curve 1) and by method II shows the results of computations using method I (r.=O.O224,, r=16z. (curve 3), z=32t. (curve 4). An increase in stepwith z=2t., 42.,8t. (curve 2), size to 7=32r. induces a marked difference in the values of Nu (curves 3 and 4). ~=16-c.and method I with Fig.2b shows the parallel results for the Nusselt number obtained by r=t* (curve 1) and by method III with 7=2t., 4t.,8~~ (curve 2). t=16r. (curve 3) and ~=32z. (curve 4). The discrepancy between the solution with ~=32r. and the solution obtained with the explicit scheme amounts to less than 6%. TO round out the investigation, we also solved the problem for large (subsonic) Fach numbers. In Fig.3 the Nusselt number Nu(t) is plotted against time for a time step-size t=0.025 and Mach numbers M=0.5 (curve I), M=0.8 (curve 2) and M=l (curve 3). At large Mach numbers the shape of the Nu(t) curves remains unchanged. Nu(t) curves correspondThe accuracy ing to larger Mach numbers lie below the Nu(t) curves for smaller Mach numbers. Nul,-, of the computations is again assessed by the discrepancy between the Nusselt numbers the condition. ANu/Nu= for steady-state solutions: at M=M,,M=M1 and M=0.5 and Nulr-t M-l - iO%. 2%is satisfied; at the discrepancy is 4%, and at M=O.8 Thus, the results of our solution of the model problem demonstrate the efficiency of treating convection problems NU using the numerical method proposed here at low Mach numbers, 1 and also the possibility of calculating flows at Fach numbers Jz z close to unity.
I r
20
10
30 t
Fig.3 REFERENCES 1.
2.
3. 4.
5.
6. 7.
8.
MAKHVILADZE G.M. and SHCHERHAK S.B ., A difference scheme for the numerical investigation of non-steady two-dimensional motions of a compressible gas. Preprint 113, Moscow, Inst. Priklad. Mekhanika Akad. Nauk SSSR, 1978. KUZNHTSOV A.E. and STRELETS M.KH., Numerical modelling of subsonic steady non-isothermal Chisl. Metody Mekh. Sploshnoi Sredy flows of a homogeneous viscous gas in channels. (Novosibirsk, vychisl. Tsentr i Inst. Teoret. Priklad. Mat. Sibirsk. Otdel. Akad. Nauk SS.SR), 14,6, 97-114, 1983. ZHMAKIN A.I. and MAKAROV YU.N., Numerical modelling of hypersonic flows of a viscous gas. Dokl. Nauk SSSR, 280, 4, 827-830, 1983. Moscow, 1961. PENENKO V.V., Methods of Fodelling Atmospheric Processes, Gidrometeorisdat, GOLOVIZNIN V.M.,KRAYZlSHKIN I.E., RYAZANOV I.A. and SAMARSKII A.A. TWO-dimensional Completely conservative difference schemes of gas dynamics with spaced velocities. Preprint 105, Moscow, Inst. Priklad. Matem. Akad. Nauk SSSR, 1983. MARCHUK G.I., Methods of Computational Mathematics, Nauka, Moscow, 1977. SAMARSKII A.A. and POPOV YU.P., Difference Methods for Solving Gas-Dynamics Problems, Nauka, MOSCOW, 1980. non-steady Navier-Stokes POLEZHAYEV V.I., Numerical solution of a system of two-dimensional equations for a compressible gas in a closed region. Izv. Akad. Nauk SSSR, Mekh. Zhidkosti i Gaza, 2, 103-111, 1967.
Translated
by D.L.
0041-5553/88 $~O.OG+O.OO 01990 Pergamon Press plc
A NUMEhICAL SIf:lULATION SCHEKE FOfiSUBSONIC FLOWS IN A VISCOUSCOMPRESSIBLE GAS*
V.A. GONCHAROV,
V.N. KRIVTSOV
and A.A. CHAFAKHCH'YAN
A difference scheme is proposed for two-dimensional flows of a compressible gas, comparable in efficiency with schemes for the equations of an incompressible fluid. 'Ihe efficiency of the scheme is determined by an original iterative process, using a procedure for eliminating for the equations of an incomthe pressure, which is characteristic A comparison with other schemes is made for a model pressible fluid. convection problem. When one is studying subsonic time-dependent flows of a gas under conditions of large temperature drops, the incompressible fluid approximation may prove unsatisfactory. Such flows are described by the full Navier-Stokes equations. It is natural to require that a difference scheme for these equations perform well even in the limiting case of practically incompressible flow - at least as well as the standard difference schemes for the equations of an incompressible fluid. This means that: 1) the volume of computational work per time step should be comparable, and 2) computations at the same time step should achieve the same The scheme proposed below meets these requirements. accuracy. In difference schemes for the incompressible fluid equations, each time step as a rule involves the solution of a difference analogue of Poisson's equation. Among the schemes available for the equations of a compressible gas are those that completely split the process with respect to the space variables (see, e.g., /l/j. Such schemes are even more efficient per time step than schemes for the equations of an incompressible fluid. However, it turns This is due to the need to calculate with a out that the second condition is not satisfied. time step-size T-h/u, where h is the grid step-size and u is the velocity of the gas. Corresponding to this step-size, however, one has large values of the Courant number (rclh, where c is the velocity of sound), which implies an increase in that part of the approximation This was shown graphically in error introduced into the scheme by the splitting procedure. /l/, where computation with a time step-size 16 times greater than that required by the condition produced a very distorted solution. Courant The difficulties involved in devising an efficient scheme for the subsonic flows of a compressible gas within the limits of the full Navier-Stokes equations have inspired the development of an approximate theory (see e.g., /2, 3/j, which permits a simpler construction This approach has undoubtedly been of difference schemes with the required properties. successful when applied to various problems. Nevertheless, we do not believe that this precludes the development of new computation schemes for the full equations for a compressible flows of a gas. Suffice it to mention that the application of the theory of time-dependent heat-conducting gas has not been adequately substantiated, as the pressure variation in such flows need not be of the same order of magnitude as the square of the Mach number.
1. We shall consider the havier-Stokes equations for a heat-conducting gas, including the effect of the force of gravity with a constant acceleration g. The system is completed by the equation of state of an ideal gas P=RpT, where P is the pressure, p is the density, The equations will be presented in a form T is the temperature and R is the gas constant. suitable for the convection problem. Let Z,pO, T0 be the characteristic length of the region, Together with the total pressure P we introduce the density and temperature, respectively. where Pst is the hydrostatic pressure of a gas of density excess pressure p=P--RpoTo-Pst, flow velocity is defined by the formula V=(gZ)"'. We introduce nonPO. The characteristic dimensional variables r&r’,
u=l.u’,
t=t’Z/V,
T=T,T’,
(P, p, P&=
poV2(P’, P',r&L Here r is the radius-vector of a point of the space, u is the velocity vector and t is the variables are indicated by primes, which will be omitted henceforth. time; non-dimensional The Reynolds, Prandtl and Mach numbers are defined by the formulae Re=p,I'Z/y,
Pr=l*C,lh,
MZ= V*/~I?T,,
y=c,/c”,
where ).Iand A are the viscosity and thermal conductivity, respectively, which will be assumed constant for simplicity,and C, and C. are the specific heats at constant pressure and constant volume. *zh.vychisl.~at.mat.Fiz.,28,12.1858-1866.1988 177
178
The transport operator in the equations can be reduced to various forms by identical For the energy equation we shall use the divergent form of the transport transformations. operator, and for the equations of motion, the skew-symmetric form (see /4/). This necessitates w=p"u. In the equations presented below the introduction of a new unknown vector-function we have used both variables u and W;this is dictated by the structure of the difference scheme 11-p'". After these transformations, We shall use the abbreviated notation proposed below. the equations are written as follows:
q $i+K(u,w) C
=-gradp+(p--l)l+ I
A(m/n)-i-$graddiv(w/n)
P-p+p,
$+div(r~w)=O, Do=-=
Irz
M*y(r-l),
C,To
1,
t+D-DpT,
D,=&,
D=+=(~M’)-I.
Here 1 is the unit vector in the direction of the gravitational force, Q, is the work of the is a skew-symmetric operator, for which we have viscous forces, K(u, w)=(ugrad)w+0.5mdivu This property of the operator K will also be carried over the identity wK(u,w)=div(uwz/2). to its difference approixmation (see /4/j, so that the schemes obtained for the equations of motion will involve a transport operator that conserves kinetic energy. The boundary conditons are as follows: for the velocity - no-slip conditions, and for the temperature - a condition of the first kind on certain sections of the boundary and thermal insulation on the others. the 2. We shall consider two-dimensional plane or axially symmetric flows, with P Differences between these two cases will be radial coordinate and s the axial coordinate. indicated in the formulae by the index v:v=O for plane flow and 11=1 for axially symmetric In these variables the differential operators are flow. grad =
alar IIa/a2 II1
divw=1dr'w,+aw, TV ar az ’
AT+-rrv$+$)T,
To fix ideas we shall assume that the force The flow region is a rectangle in the r, z plane. In the I‘, s plane we introduce a rectangular grid of gravity acts along the s coordinate. with constant step-sizes h, and h,, whose points are numered by integer indices i, j. Following /5/, we define the thermodynamic variables at the centres of the grid cells:(P, T, P)<+~,i+~h~ and the components of the vector W at the midpoints of the appropriate sides: WV,,. i+'br zaz, ‘+‘h. J. The values of the difference functions at other points of the grid are obtained by Linear interpolation from the nearest grid-points at which these functions are defined. The boundary of the region consists of several grid lines, on which the normal velocity To specify the temperature component is defined, as well as the thermal flux, if necessary. and the tangential component of the velocity vector on the boundary, we introduce fictitious cells, whose data are defined in such a way that linear interpolation together with the data in neighbouring cells gives the prescribed boundary value. In constructing the difference scheme, we shall use the idea of splitting wherever we consider it justified. To obtain second-order time accuracy, we will use the symmetric splitting scheme of /6/. The sequence of computations per time step is determined by the following difference formulae: (w,),+0.5K,h(U,, wr"")=O, (pT),f0.5
(w,),+0.5K,"(u,, f1'*05)=o,
div,"[u(pT)'.']=O,
(w,),+O.5K,"(u,, w,o")=O,
Re-’
(PT),=D,A,“T’,~,
A,‘“( w,O “/a),
(p),=O;
(lb)
(w,),i-0.5K.h(&, w,"~)=O,
(pT),+0.5div,"[z~(pT)~~j=O, (tc,),=O.5q-’
(p),=O;
(P) r=O; (w,),=O.5q-’
(la)
(2a) (P)
Re-’ A,*( to,” “/q).
(3rd (3b)
179
,
(2~),=0.5q-~Re-'A," (pT),=D,A,"T",
(~,).=0.511-'Re-'A,L(~),
q
(p).=O;
(4b)
(w,),=-$~rad,"div'(~),
(,?a)
( 1 1;
grad,"div"E
(wz)T=&
(5b)
$=[p(pO.J,T'.")-l]Z,,
n(w),=-grad"p"+$,
(44
(pT)~=--D,~0.5diP(~a5/7])f~yh,
(p)z+div"(qw)=O,
(6a) (6b)
P=p+P&.&D=DpT. Computations using formulae (l)-(5) are then carried out in the reverse order. Here (v)z=(qn(P~)/T, where T is the time step-size, 'pn is the "new" value of the difference function, to be determined at the current step of the splitting scheme, and 'p. is the "old" value, known from the previous step; @."=0.5(cpnt~o). The function p(p, T) at step (6a) is determined from the equation of state. The difference functions u and n are computed from the data on the previous time layer and not changed on the current time layer. The superscript h denotes difference differentiation operators with respect to the appropriate variable, the latter being indicated by subscripts r or z. These operators have a natural form for the grid: Mlh(&f)
[&‘“fl
lIJ+‘b=
[ Vu) *+‘/,fi+*(Jw I-*k-1 lj+1,, I Zh,r,‘l
iW)r-(rvl)r-l ,,,+'I>= -h,= I I-‘ ' 3+ !- L, 1,+'I* 1
(Pf)S+l-(rVl)
[rdl (fn++-f,+~IJ) -riy (f,+~,,-f,_s,JI, [A:fl z+a/>,, = r:+,,xh.* (A.“f)i+v,,,=
f+%[f~+I-2fj+fj-i1 h.’
’
f[fj+a,:-2f,+a+fi-al (A.“f)t.,+‘,:= It,* ’
km&” /Ir+‘,>.j =
fi+‘b.j+lh-f,+‘l.,l-‘h~ hz ’
f f+B,~+l-ff+‘lr., (div,”f) f+‘hJ+% = 11, ’ div” w* = div,"LU,+ div,"w,. ( I(', 1 In some of these formulae subscripts common to several variables have been written outside the brackets (the first subscript before the brackets, the second after them). The expression for the work of the viscous forces also has a natural form, not shown here (in convection problems this term is small compared with the others). The approximation of the operators K(tr,U?) was constructed as in /4/. Steps (l)-(4) are implemented using one-dimensional pivotal condensations. Step (5) is accomplished by a version of the method of alternating directions. The function wI in Eq. (5a) is assumed to be known. Then ZL',is found by one-dimensional pivotal condensations and substituted into Eq.(Sb) , which is solved by one-dimensional pivotal condensations for mz, and so on. for In practice it has been found that two or three iterations are sufficient convergence. We will now consider the system of algebraic equations arising at step (6). The unknowns 10,p, T,p. We use a well-known device /7/, in which the here are the difference functions equations are divided into two groups: the first includes the energy equation, and the second, the remaining equations. First, assuming T to be known, one solves the equations of the
180 second group and finds the other unknown functions. Then T is determined from the energy equation, and so on. In performing these computations it was assumed that the maximumadmissrelative change in temperature in one extraneous iteration was ible 6=1o-s; to achieve this it wassufficient to perform one or two extraneous iterations. The efficiency of the scheme is largely determined by the iterative process for solving the second group of equations, which will be consideredinSect.3. 3. We introduce a new difference vector function V'T)IlF. Then the system of equations of step (6), given the difference function T, is written in the form 2(v),=-gradhpO,'++,
(uL=(v-_gWo)l',
(7)
(P)~=-div"v,
(3)
p=DpT--D---Pst
(9)
The function $ in the equation of motion will be considered to be known. We first eliminate the difference functionspand p. To this end, we substitute (9) and substitute the result into (7). This gives 2(~),=0.5D~ grad"T div” v+q, (P-Q-0.5
gradh(po--Ps-,-0,5D
(8) into
(1Oa) gradh(po7’),
(lob)
where the subscript "o", as before, denotes functions known from previous steps of the splitting scheme. Rather than solve this system directly, we shall eliminate the pLessure from the equations of motion, a normal procedure when solving the equations of an incompressible fluid. We introduce the difference operator rothu=grad,'hu,-grad,'hu,, where grad."'and grad,'* are the operators of difference differentiation with respect to the respective variables, whose values are defined at the integer-indexed points of the grid: (grad,.hf),,
=
fy-fw,
(grad,.“f)
,j
=
h+;-h--y
I
It is readily Applying
seen that rothgradh=O. the operator rot*
to Eq.(7),
we obtain
rot"u=f,=lT rothIp+roth u. (11) Assuming now that (P)~ in (8) is given, we obtain a difference analogue, of the well-known problem of determining a vector field, given its divergence and curl. If the normal component of the vector on the boundary of the region is given, theproblem is not necessarily solvable for any function defining the divergence: the integral of the function must equal the integral of the normal component of the vector along the boundary. In the case of the difference equations (8) and (11) this property of the differential problem is represented by the fact that the number of equations (for a given function (p)* exceeds the number of unknowns by one. If the function (P)~ is properly chosen, the rank of the system of equations (8), (11) equals the number of unknowns, since the sum of all the equations of (8) becomes an identity, inthis case Zdi@u=O, where the summation extends over all cell centres and it is taken into consideration that the normal component of the velocity vanishes on the boundary. Below, (P)~) will always be caiculated from Eq.(8) for some given function V satisfying the boundary conditions, and therefore the system of Eqs.(8), (11) is consistent. We will now construct a procedure for eliminating one of the components of vfrom systen(8), (11) making use of the analogy with the familiar formula rotrot=graddiv-A for the differential operators. We formally extend the definition of Eq.(ll) to the boundary lines of the grid r=conSt, say by setting f,=O. To this end it is necessary to introduce fictitious cells near the boundaries and suitably define the function u, in them. Since v,=O on the boundary, we arrive at the condition VI,%,j=JJ., ‘i:, ,, (12) where the subscript '/*corresponds to the fictitious cells and '/z corresponds to the nearest interior cells. We will now consider the operator div?"', whose values are defined at the same points as the component v,,with the exception of the boundary lines z=const: (div,'" f) 1+%.,=
ff+l,jl~:+L--fUTiV “:+;, hr
.
It is readily seen that div,"'gradr*"=gradZhdivlh. Apply the operator operator grad," to ~q.(8) and then subtract one from the other. Poisson equation for v,: A."u,=f=-gradZh(p),-div,'"f,, A"=grad,"div,"+div,'h grad,'h.
div,'" to Eq.(ll) and the This gives a "difference" (13)
181 and a condition of type (12) on the lines The boundary conditions are U.-O equation was solved using the direct method of the ordinary Fourier transform
z=const. This /6/.
NU
4
b
a t
: 3
2 : 31 k 0
10
20 Fig.1
NU
4 t
a
Fig.2 We can now formulate the proposed iterative process for solving Eqs.(i')-(9). Fix some Substitute the resulting function function (pf=, say identically zero, and solve Eq.(13). uI into (loa), and then solve this equation for U, using one-dimensional pivotal condensations. pivotal Kow substitute the resulting function V, into (lob), and again use one-dimensional condensations to find a new v,. The last step of the iteration is to substitute the computed vector U into Eq.(8) and calculate a new function (p)%. In all our computations, the iterative process converged fairly rapidly - in two or three iterations - even for comparatively large Mach numbers (0.8 to 1.0). Counting also the extraneous iterations due to the change in temperature, one time step required solving Poisson's equation at most six times. As shown by comparison, the computation of strongly compressed flows takes all in all 2 to 2.5 times m3re time than the analogous computation with the Boussinesq approximation. 4. Th e proposed numerical method (method III) was compared with the explicit scheme of /6/ (method I) and the implicit splitting scheme of /l/ (method II) for the model problem of /l/. Consider the motion of a compressible gas in a closed square of unit length under the action of the force of gravity. The upper and lower walls are thermally insulated, but a constant temperature is maintained at r=l , T= 1.5 at on the side walls: T=i r=O. The velocity satisfies the no-slip condition at the boundary of the domain: u,=u,=o. The initial conditions are as follows: a linear temperature distribution in the gas at rest: T=l.S-0.5, uI,u,=O, p=D=const, corresponding to the steady state in the absence of gravity. bon-dimensional parameters: Re=250. Pr=O.71, y=1.4. The problem is solved for two values of the Mach number: M,z=1/320, MZ2=1/20. The result of the computation is the behaviour with time of the mean Nusselt number at the hot wall:
The computations were carred out with different time step-sizes on a 21x21 grid. To check the accuracy of the computations in the steady section, it was verified that Nuj,,o= NIlI,-,. In Fig.la the Kusselt number Nu(t) at the hot wall is plotted against time for M=M,, 1.=0.0056 and method II with r=2r., 42.,8% using method I with maximum admissible step-size (curve 2) and ~=32~. (curve 3). An eightfold increase in step-size (curve 11, r=16r. yields a deviation of less than 1.5% from the solution of the explicit scheme, and therefore A further increase in all these results are represented by the same curve on the graph.
182 step-size caused a marked change in the Nu(t) values (curves 2 and 3). Fig.lb shows similar olots of the Nusselt number Nu(t) obtained by method I with r=r.(curve (curve 2) , r=32r. (curve 3), z=64r. (curve 4) 1) and by method III with r=&., 4t., 8r., 16~. and r=128~. The solutions represented by curves 1, 2 and 3 differ by less than (curve 5). 2%. In computations with .M=M, using method III, the discrepancy between the mean Nusselt numbers of the cold and hot walls amounts to ANu/Nu=2%. This discrepancy determines the accuracy of the computations. Thus, the solutions corresponding to curves 1, 2 and 3 are identical, to within the accuracy of the computations. The discrepancy between the solution to -5%. The larger with ~=64r. and the solution obtained by the explicit scheme amunts difference between these solutions is obvious1 y due to the fact that the steady solutions of the splitting scheme depend on the time step-size 'c. It is also worth noting that all the curves obtained by method III duplicate the shape of curve 1, whereas the curves produced by method II show considerable fluctuations about the solution obtained from the explicit scheme. Similar results were obtained for the other value of the Mach number, M=Mt. Fig.2a curve 1) and by method II shows the results of computations using method I (r.=O.O224,, r=16z. (curve 3), z=32t. (curve 4). An increase in stepwith z=2t., 42.,8t. (curve 2), size to 7=32r. induces a marked difference in the values of Nu (curves 3 and 4). ~=16-c.and method I with Fig.2b shows the parallel results for the Nusselt number obtained by r=t* (curve 1) and by method III with 7=2t., 4t.,8~~ (curve 2). t=16r. (curve 3) and ~=32z. (curve 4). The discrepancy between the solution with ~=32r. and the solution obtained with the explicit scheme amounts to less than 6%. TO round out the investigation, we also solved the problem for large (subsonic) Fach numbers. In Fig.3 the Nusselt number Nu(t) is plotted against time for a time step-size t=0.025 and Mach numbers M=0.5 (curve I), M=0.8 (curve 2) and M=l (curve 3). At large Mach numbers the shape of the Nu(t) curves remains unchanged. Nu(t) curves correspondThe accuracy ing to larger Mach numbers lie below the Nu(t) curves for smaller Mach numbers. Nul,-, of the computations is again assessed by the discrepancy between the Nusselt numbers the condition. ANu/Nu= for steady-state solutions: at M=M,,M=M1 and M=0.5 and Nulr-t M-l - iO%. 2%is satisfied; at the discrepancy is 4%, and at M=O.8 Thus, the results of our solution of the model problem demonstrate the efficiency of treating convection problems NU using the numerical method proposed here at low Mach numbers, 1 and also the possibility of calculating flows at Fach numbers Jz z close to unity.
I r
20
10
30 t
Fig.3 REFERENCES 1.
2.
3. 4.
5.
6. 7.
8.
MAKHVILADZE G.M. and SHCHERHAK S.B ., A difference scheme for the numerical investigation of non-steady two-dimensional motions of a compressible gas. Preprint 113, Moscow, Inst. Priklad. Mekhanika Akad. Nauk SSSR, 1978. KUZNHTSOV A.E. and STRELETS M.KH., Numerical modelling of subsonic steady non-isothermal Chisl. Metody Mekh. Sploshnoi Sredy flows of a homogeneous viscous gas in channels. (Novosibirsk, vychisl. Tsentr i Inst. Teoret. Priklad. Mat. Sibirsk. Otdel. Akad. Nauk SS.SR), 14,6, 97-114, 1983. ZHMAKIN A.I. and MAKAROV YU.N., Numerical modelling of hypersonic flows of a viscous gas. Dokl. Nauk SSSR, 280, 4, 827-830, 1983. Moscow, 1961. PENENKO V.V., Methods of Fodelling Atmospheric Processes, Gidrometeorisdat, GOLOVIZNIN V.M.,KRAYZlSHKIN I.E., RYAZANOV I.A. and SAMARSKII A.A. TWO-dimensional Completely conservative difference schemes of gas dynamics with spaced velocities. Preprint 105, Moscow, Inst. Priklad. Matem. Akad. Nauk SSSR, 1983. MARCHUK G.I., Methods of Computational Mathematics, Nauka, Moscow, 1977. SAMARSKII A.A. and POPOV YU.P., Difference Methods for Solving Gas-Dynamics Problems, Nauka, MOSCOW, 1980. non-steady Navier-Stokes POLEZHAYEV V.I., Numerical solution of a system of two-dimensional equations for a compressible gas in a closed region. Izv. Akad. Nauk SSSR, Mekh. Zhidkosti i Gaza, 2, 103-111, 1967.
Translated
by D.L.