- Email: [email protected]
A method for the numerical investigation of unsteady subsonic flows of viscous gas in ducts
USSR Comput. Maths. Math. Phys. VoL 21, No. 5, pp. 143-159,198l. Printed in Great Britain
0041-5553/81/050143-17$07.50/O 0 1982. Pergamon Press Ltd.
A METHOD FOR THE NUMERICAL INVESTIGATION OF UNSTEADY SUBSONIC FLOWS OF VISCOUS GAS IN DUCTS* A. T. FEDORCHENKO Moscow (Received 14 February 1980)
A NUMBER of computational models for solving non-linear problems of unsteady subsonic flows of a viscous gas in a finite segment of a duct are presented. A difference scheme is presented for the integration of the complete Navier-Stokes system of equations for a compressible medium. Some general features of the application of mesh methods for investigating three-dimensional acoustic processes in ducts are considered. Examples are given of the realization of computational models and of the difference scheme in a series of unsteady problems of the development of the axisymmetric flow of a viscous gas in a short cylindrical duct.
Introduction At the present time in the general development of computational methods of gas dynamics the method of computing essentially unsteady subsonic flows of a viscous gas in ducts using the multidimensional Navier-Stokes equations is still insufficiently developed. The absence of solutions of non-linear boundary value problems of this class is objectively explained both by the difficulties in constructing models of unsteady subsonic flows of a viscous gas in bounded regions (this is due to specific requirements imposed on the boundary conditions), and also by problems of the correct resolution of the three-dimensional acoustic processes by a difference scheme. In particular, the investigation of such topics would facilitate the solution of important applied problems in non-linear acoustics [ I 1. In this paper a group of models is presented for computing the unsteady subsonic flows of a viscous heat-conducting gas (using the complete three-dimensional Navier-Stokes system of equations) in a finite segment of an arbitrary duct. The acoustic characteristics of penetrable end barriers of various types are investigated, which enables methods of efficient action on the longitudinal oscillations in the duct to be indicated (for example, for the purpose of damping the low-frequeny components of the spectrum). For application to a wide class of gas-dynamic problems a difference scheme for integrating the complete Navier-Stokes system of equations for a compressible medium is proposed. The fundamental characteristics of the scheme are given. Some necessary mesh criteria for the solvability of the acoustic processes in the ducts are determined. As examples of the use of the proposed models and the difference scheme numerical solutions are presented of a series of unsteady axisymmetric problems in which non-linear oscillatory processes in the development of the flow of an originally quiescent viscous gas in a short cylindrical duct are considered. *Zh. vjkhisl. Mat. mat. Fiz., 21,5, 1215-1232,
1981.
143
A. T. Fedorchenko
144
1. General formulation
of the problem
Consider a bounded domain G, which is a finite segment of a duct with rigid impermeable side walls. For the general formulation of the problem we introduce some curvilinear orthogonal coordinate system {OX. Oy, 02). assuming that the axis Ox is directed along a slightly curved nominal axis of the duct. We assume for simplicity that the permeable end boundaries rj, j = 1,2, are elements of the coordinate planes x = xl, x =x2, x2 - x1 = L. The duct wall, the boundary r,- (I’==~,uI’,-ur2, I’,,,=r,, n I-,,, C=GL r) , is defined by some continuous piecewise-smooth surface _!Z(2, zj, z) =O. The subsonic unsteady flows of a viscous heat-conducting gas are modelled numerically in the domain G. Without significant loss of generality we wilJ use the model of a one-component ideal gas, governed by the equation of state p = pT/-y, y = cPo/c,o = const, which is in equilibrium in the conditions of the low-frequency approximation. Then, for the five independent parameters of the flow ~1, ~2, ~3, p, T the initial system of five differential equations includes the continuity equation, three equations of motion (the complete Navier-Stokes equations) and the energy equation. These equations are complemented by the local relations p = p (T, p), X=X (T, p), . . . . Here and below we consider the dimensionless values of the coordinates (x, y, z), the velocities (~1, ~42,ug), the time r, the density p, the pressure p, the temperature T, the coefficients of shear viscosity p (the coefficient of volume viscosity is ignored) and the thermal conductivity x. These quantities are expressed in the units rg, co, r&g, po, poco2, To, ~0, x0 , respectively. The characteristic values chosen were rg the characteristic half-width (radius) of the duct, To the characteristic temperature of the walls, co (To) the adiabatic velocity of sound, po, ~0, x0, cP,, , the density, viscosity, thermal conductivity and specific heat at the temperature To and some characteristic (initial) pressure in the duct. We define the fundamental similarity criteria:
In the domain G for the interval R = (0, tk) we can formulate the unsteady boundary value (more accurately, mixed) problem with the initial field {u, p, 2’)~ at t = 0 and boundary conditions on r. On the impermeable wall of the duct rw, conditions of the following form are usual:
u=o,
a,dTldn-I-a?T=a,, ai=ai(r, 1).
d,,
t4.
(1-l)
Some element of the boundary may be a plane (axis) of symmetry, where the conditions
u,,=o,
du,ldn=dpldn=dTldn=O,
(1.2)
are specified, and u,, U, are respectively the normal and tangential components of the velocity vector. The specific models of the flow in the domain G discussed below are primarily defined by the method of specifying the boundary conditions on rl, r2.
Unsteady subsonic flows of viscousgas
145
2. Basic computational models The proposed non-linear models of unsteady subsonic flows in the domain G include as important elements two types of permeable boundary surfaces (we denote them by S,,,*, S,*), intended for specifying the conditions on ri. To be specific we assume that at all the points r E T’i considered we have U, = (u, n) # 0. The upper suffvres “t”, “-” denote the intake of the flow at the given point T’i: (u, n) > 0, and the outflow: (u, n) < 0, n is the inward normal. 1. Let us first consider the &-surface (some simplified versions of this model were proposed previously in [2,3]). By the general model of S,, at each point r E S, + C Pi, t E 52, conditions are used of the form
\‘(r. f)=f(r.
t, p(r. t)).
7’!r, t) =g(r.
f-{:i+
(2.1)
!:> in}.
t, p(r, t)),
(2.2)
where f, g are given functions. Here and below Vi = pF ui (6 is a specified real parameter, in general different for i = 1, 2, 3). Sometimes instead of (2.1) it is preferable to specify the function $ and the unit vector b: P’luI =$(r,
L p(r, t)),
u=lulb(r,
t).
(2.3)
In the case U, < 0 (on S,,, - ) a single condition is imposed of the form
t’Tl(r, t)=(V.
n)=f),(r,
t, p(r,
t)),
rES,, -=I?,,
J&.
(2.4)
We will discuss separately conditions of the type (2.2). In the simplest case, ag/ap = 0, the temperature distribution
I’r=g (r, t) ,
rEr,+,
t4.
(2.5)
can be specified on ri+. Physically this condition corresponds not to the cross section of the developed flow, but rather to a permeable wall (lattice), and therefore it would be more suitable for a permeable segment of the boundary T’,+ Obviously, even for the whole boundary Pi = Pi+ it is possible to use the more general relation T(r, t)=gl(r,
L) [P(r,t)]‘:-‘“‘+g?(r,
t),
glao,
(2.6)
which may be equivalent to the specification of a slowly varying entropy profile s (r, t) over a large part of the cross section Pi+ (except for a neighbourhood close to Pwi with characteristic width of the order of the thickness of the acoustic boundary layer 6, = max (SF’ ,a%’ > along the wall of the duct, see (4.7)). We note that in a number of cases the use of (2.5) instead of (2.6) should not lead to significant differences in the solutions, at least for the low-frequency band of oscillations (see below remark 2 and subsection 3 in section .5)., We present some particular cases of the model S,.
146
A. T. Fedorchenko
Within the limits of the S ,+,-model the instantaneous values of V, (r, t) at the point r E F’j are determined (with possible time delays ri) by the drop between the pressure p (on the “mternal” side of I’i) and the “external”
pressure pe:
v,, (r, t) = cp(r, t) I AP Ik sign (APL AP =
:i
s:,
me (r, tl) -
f2 (r, t, VI (r, t)),
vt (r, t) =
T (r, t) = g (r, t,
i
where r,>O,
920,
p,>O,
0, t),
i=l,
2,
V3 = fs (r, 4 VI (r, th
(2.7)
P Cry W;
t) 1bplk sign (Ap),
V,,(r, t) =q(r,
Sm,-:
k > i!i=t-Ti(r,
P (1, tz),
g>O, and
In some cases the distributions
f2, f3 are specified functions.
p (r, t), pe (r, t) can be replaced by the average values po (t),
p,o (t) over the cross section ri:
v, (r, t) =(r
(r, 0 I APO Iksign
(Ap”),
ApO=p.O(t,)
... .
-p”(t,),
(2.8)
In the simplest version of the model S, (we denote it provisionally by S,,, 0) it is assumed that the vector V at the given point r E rl .+ is specified, that is, it is independent of the oscillations of the pressure inside the domain G: Smo+: V (r, I) =f (r, t),
T(r,
t),
Smo-: V,(r, t) =fn(r,
t)=gk,
6 p(r, t)), (2.9)
tEQ2.
rEl?j,
The model Snto+ can be formally regarded as the result of the passage to the limit in Sm+,+: pe + 00, cp+ 0, p < const. However for S, 0 - such an interpretation
is invalid (pe 2 p for (u, n) < 0).
2. By the model S,, for r E Pi, r E R there may be conditions
T(r,t)=g(r,t,p(r,t)),
~(r,t)=f,(r,t,V,(r,t)),
S,+. -1
&-:
of the form
u/lul=b(r,t): p (r, t)
(2.10)
=fP b, 4 Vn(r, t) ) ;
fp, b, g are specified functions. In a particular case (the model S,n) af,/aV, = 0. In &other
particular
S,+: p(r, SP-:
p(r,
version of the Sp-model the conditions
t)=fr(t, l) =fP(t,
mj(tj)), W(tj)
T(r, )7
t)=g(r,
r=rf,
4 p(r,
t)),
b=b(r,t),
(2.11)
t42,
are possible, where fi = t - rj (t), ~j > 0 is a specified delay parameter,
rni is the instantaneous mass outflow in the cross section I’i (in calculations for ri = 0 the approximate values of rnjdetermined in the cross sections x1 + hX or x2 - hX, 0 < hx < 1 may be used).
In some cases (for example, if a~p+azV’,=a,, equivalent
to S,.
ai=ar(r,
t) , rd?,)
the model SP may be
Unsteady subsonic flows of viscousgas
147
The models S,, S, can be naturally generalized and made more complicated by introducing “controlling boundary conditions” of the type
V(r, t)=f(r,
t
)...,
Jr),
rdm+crj,
Ed&
(2.12)
where W’ is a finite set of parameters of the process investigated from the solution of the problem on the point set G’XQ’, G’cG, R’c [t,, t) d2. Similar analogs are possible for fP, fn,g,b, . .. also. In particular, conditions (2.7), (2.11) with ri > 0 can be reduced to the simple form (2.12). In general, the process of numerical modelling of the gas-dynamic phenomena in the domain C inescapably involves the need for a continuous investigation in the computational process of some set of parameters of the process of solving the problem in the preceding time interval (with local estimates of the number M and the sign of (u, n) in I-j, with an analysis of the spectrum of the acoustic oscillations generated and selection of the next algorithm of operation on the oscillatory process etc.). Remark 1. The assumption of the subsonic nature of the flow in c (in every case, M < 1 close to I’) imposes certain constraints on the functions f, fp, g, . . . in (2.1) - (2.11). In some simple cases (see section 5) preliminary estimates of the range of parameters in the solution of the problem are possible. If it is impossible to provide such preliminary estimates the tracking of M and also of other parameters in G may be accomplished by using controlling conditions of the type (2.12). 4.
An approximate analysis of the initial system of differential equations close to rj enables the number of independent boundary conditions in the models S, +, SP* to be validated. Indeed, let the boundary Fj be a normal cross section of a quasi-parallel subsonic flow. Then near each point rj the initial system of five equations can be interpreted as locally hyperbolic in the variables x, t (see, for example, [4]). Here, all the gradients along the axis Ox are assumed to be small, and the terms with derivatives along Oy, Oz are interpreted as distributed sources or force functions. A subsequent characteristic analysis [4, 51 of the system leads to the result that the required number of boundary conditions at each point r E I’j+ must be N+ = 4, and at the point r E rj- it will be N- = 1 (correspondingly, for the two-dimensional problem in x, y coordinates we have N+ = 3, N- = 1). It must be stipulated that a similar result is invalid for estimating the number of boundary conditions on a rigid permeable wall of the duct rw. In this case with the dominating role of the parabolic processes for ~2, ~43, T (with respect to x, f in the equations of motion and energy) the number of boundary conditions may amount to N+ = N- = 4, although here also the basic elements of the models S,, SP can be used to determine the connection between ul, p, T. It is important that the very nature of the proposed system of boundary conditions (2.1) (2.11) differs significantly from the traditional approaches of various authors (for example, in [4] for I?,+ it is proposed to specify the angles of approach of the flow, the stagnation temperature and the stagnation pressure, that is, conditions obviously inadmissible for the modelling of viscous subsonic flows, especially for taking into account non-linear acoustic processes).
148
A. T. Fedorchenko
5. It is natural that in (2.1) - (2.11) the number of independent boundary conditions is always less than five. Therefore, in the numerical solution of the problem the parameters at the boundary points must be calculated by introducing additional difference algorithms, containing to some degree elements of spatial extrapolation from the internal points of the domain. These algorithms are based, in particular, on the characteristic dependence of the pressure at the points r E S, C r. I E f’,+, on the process in the entire domain G by means of the acoustic connection ? (M < 1 in G). For example, when the difference scheme of Section 4 is used the pressure ,;+r (r=S,, r=r,, tn+‘= (n+l) ‘c) can be computed by normal extrapolation of second-order accuracy from internal points of the domain G, already computed on the layer tn+l (sometimes this procedure is carried out by using a difference analog of the continuity equation, written in conservative form for the boundary cells of the mesh). Then conditions (l.l), (2.1), (2.2) enable us to compute VF+l and T{” on S,+, &?+l, $?+I,. . . . Various extrapolation procedures (for example of the type (d2~!dn2)ri+i=O)can be used to compute Tf!‘l on S,,- and S,-. V;+l on S,+. 6. We will determine the general model of subsonic flow in the domain G as a combination of the S, and SP surfaces used to specify the boundary conditions on l?j. Therefore, three basic types of computing models can be introduced: {m, p}={S&& {m,
S&)1,
(2.13)
~}={~n(~i)~Sm(~A))~
{P,P}={S,(5i),Sp(5k)},
ifk,
i=l,2,
F;=l,2.
We mention briefly some of their features and limitations. If in terms of the model (2.13) the surfaces rj are identified with the end cross sections of a complete length of duct, then in this case it is difficult to take into account local end effects close to I’,+,i(the vortex flow past angles and sharp edges, diffraction wave phenomena etc.), depending on the shape of the boundaries and the conditions of the flow outside the domain G. However, this problem can be solved if, using the same models of S,, S,, we consider the system of “ a duct inside a duct” (that is, a comparatively narrow duct connects domains with a sharp expansion of the transverse cross section). Moreover, the models of (2.13) are correct only for wave processes with weak discontinuities, since the conditions of reflection of shock waves from the boundaries I’i are more complex in nature, going beyond the limits of the models S,,,, S,. Using the model SPo it is practically impossible to specify in the “input” cross section rj (in the sense of the direction of the averaged mass flow) arbitrary profdes Vl with characteristic values of the transverse gradients (especially with a definite value of ml (r)). Therefore, in the majority of cases, except for a number of specific problems (see subsection 6 in Section 5, the model $0 can be usefully used only in the “output” cross section.
3. The acoustic characteristics of the permeable boundaries By linearizing the boundary relations (2.1) - (2.11) it is possible to estimate the acoustic properties of the permeable boundaries r’.
Unsteadysubsonic flows of viscousgas
149
1. Let us consider a small element dTj=rj (to be specific, j = l), on which, using the model (2.9), we put V,=pFui=f (t) , where 6 = const, f(t) is a function, which is fairly smooth and monotonic for t E (r 1, t2)c Cl, and is practically unchanged within a characteristic period of oscillation fp = 2rr/w < t2 - rl . Close to dl?l* (on clTl+ for g2 = 0, gl G=const > 0 in (2.6)) we will consider small normal adiabatic oscillations (u,=u=ii+u’, p=pi-p’, p =pSp’, p’=p’/~‘, / u’iii 1
&i/C=---jii?
sign (li), (3.1)
‘I= (l+q)/(l-q)
for q*1.
It is characteristic that 9, n are real and are independent of the frequency of oscihation. We note that 9 ~-l(~~O)f~r~~~/ii(thatis,1~1>1i;9-+O(~~1)asC;~O.For~E(O,~~/~~)we have 9 < 0,~ < 1 on d rl+ and 9 > 0, n > 1 on d I’1 -. The case j = 2 can be discussed similarly. However, significantly different results are possible when the oscillating functionsf (t) are specified on&O C q.
2. Let us consider the normal harmonic oscillations (U‘=uo’exp (--id), p’==po’ exp (--Lot) ) close to the element d I’l, on which the model (2.7) is used with k = 1, ~1 = 0. We assume that cp(t), pe (t), ~2 (t) are practically unchanged in the course of the characteristic period of oscillation rp. Then for drl*
Putting 11~ il5.
T:=& p’Z”(c. (,LI~-~). we can rewrite (3.2) in the form
E---
q
u 2L.t c [ Pe-P
I
Fz-
Wra+j(pe-_P)) Ipp-jd
(3.3)
*
For example, 9 * 0,~ = 1 for zx-rP/(p&,-I_i). For [= 1 (VI = pul) we have Yp+c (PC-p) =p (y-1) +p20, that is, 9 < 0 on drl*. We note that (3.3) reduces to (3.1) only for the case when (u, n) > 0 as (Y+ 0, pe + m, P < const. Remark 2. Just at the points r E lYlwith the condition (2.5) for ag / dr=Owe have p’=:yp’ / T?. The relation p’=p’ / ~2 is true only for rs.r,t 6:’ (see (4.7)), which by the way does not change the order of the estimates (3.1) - (3.3). In any case the portion of the acoustic energy absorbed in such a case due to heat exchange (see [6,7]), can be estimated by the quantity AIU/ W= 0 (0, / R&Pr)‘+ usually extremely small in the low-frequency band.
The total effect of the weakly curved surface J’i (with the distributed quantity 9 (r, t), I E T’i)can be estimated approximately by “averaging over the cross section” the values of the conductance Qj and the reflection coefficient nj:
Qj(t> m
+JJq(r,t)da-(Qj),,+i(Qj)~~, ’ r1
j=1,2,
A. T. Fedorchenko
150
= 0, then practically total transmission of the quasi-plane normal wave through the If CQj>im cross section rj (that is, I qj I < 1) must be observed if (Ql Xe = -1 or (Q2)re = 1. Putting (Ql),, > 0 or (Q2)re < 0 may lead to the amplification of the longitudinal (along Ox) oscillations in the domain C. 3. Let us now consider small normal oscillations close to the boundary r2, on which by (2.11) we have put Pb2,
Y, 2, t> =fp
0) i-P@) i.m (0 -6
t’=t--tm,
(t’> I,
Sm),O,
(3.4)
where
The quantity E2 can sometimes be estimated as follows: if on rl the model Smo+ with ~=l,&n~/at~Oisused,thenE~=:m~,andifonI’~ themodel,Spisused,thaniii2~(ml +mz) 12, and so on. Let us assume that fP, 13are slowly varying specified functions (that is, practically constant within a period), and the gradients aTjay, aT/az are small close to F2. Then if rm = 0 then for the quantities uO’, pO’, pO’ averaged over the cross section r2 (with area 02) it is easy to obtain p’=@?[ u,‘p,+p,‘U,]. Putting pG’~p0’,k2, B=?cT?, we have p,‘c, (I-B&l&*)
=Buo’Po&
(3.5)
Q+ (if,--Bii,/&) lB= [ &-B.U, sign (U,) l/B.
(3.6)
We note thst Q2 = 0 for B~:Co2/ii,, U,+O. For C, = 1, H, < 1 the quantity Q2 = 1 (that is, 1~2I Q 1) is attained for B = 1. By specifying the parameter r m > 0 in (3.4) we can obtain various complex values of Q2 (t, w).
4. The difference scheme 1. For the numerical integration of the initial system of differential equations we can use an explicit difference scheme of the predictor-corrector class with the property of over-all approximation, which is a generalization of the scheme previously proposed in [ 81, This scheme, satisfying a number of specific requirements, was successfully realized in an entire series of essentially unsteady problems of the dynamics of a viscous gas [Z, 3,9]. We note that, in general, the models of Section 2 are independent of the type of numerical scheme used. We briefly present the main features of the difference algorithm by the example of the model equation for the vector function w (x, t), x = (xl, x2) E G, t E 52: (4.1)
fR=f(w),
P.=P (w> f
~=O(l/Re)
k, l=l,
2.
Unsteady subsonic flows of
gas
W’SCOUS
151
We will consider some unsteady boundary value problem in the domain G with boundary I’ consisting of rectilinear segments parallel to the coordinate axes. We will construct a non-uniform rectangular mesh (formed from straight lines parallel to the coordinate axes) and consider the set Ch of nodes (z,‘, x2j) EC, eh=GhU l?,,. We will also introduce the set Gh * = Cl * U G2 * of “complementary” nodes situated at the middles of the segments of the coordinate straight lines connecting the nodes of the basic mesh: Gi’= {iPhJ} = {xi=
(~i’+x:*~
Gz*= {s.+*“‘} = {zi=z:,
) /2,
IC~=Z~~, (xi, s,) EG},
zz= (s:+zf’
) /2,
(CC,,x,) EG} .
We define the time mesh: &G
p+o= (n+a)T;
{&O,
cJ”E, 1; n=O, 1,2,.
We will denote the mesh functions at the nodes ChXor’,
..},
&E[0.5, 11.
GhXure by w, w* respectively.
Then the difference scheme for (4.1) is realized within an irregular 9-point template of ~??hin two consecutive stages: 1) at the predictor stage W;,~~~~j=(Wij”+W~*,,j)/2--EZ[DI:)’,,+(~~2’
(wjn+W:j*i)12
Wr.L=
+Di’,“:
/2+0;3f”,
--ET[ (is,“‘+S,‘::) Dj’:I,,f”=*
DI:)l,,fn=~(f,=,,j-f~jn)lhl*‘,
~~‘)fb(p
+ai+'(ftF+!,j-fij")
(f
~i(“fn=az-‘(fu”-f*;-l)
)/Z]f”,
lh,*‘,
(fC,j+l-fzl”)
(4.2)
7
+az+‘(fi;+i-fij”).
h,*’ = I2,‘~-3i
tt1
hl*L Id-& kl I,
1,
ah*‘=hk=‘/hk” (h,+‘+h,-‘) ;
2) at the corrector stage
n+1
Wij
=wij “-t[
Dd“f=2a,-’
(jT*-f&,,)
Do(2)f=bz-i(fz-fi,j-~/a>
fi=Vfijn+
=[O,
+2a,+i(fi;~/z,j-~~)1 +2~+‘(fi~j+~/z-_T*)
(3-V)
fz=Vfijn+(l-V)
(D~‘)+D0(‘))f-~gij”],
)
(4.3)
(b2-‘f~,j-l,*+b2+‘fi~]+1!I)1
(bi-‘f:-‘h,j+bi+‘fi+lil,j),
II,
t&*l=h=‘/ (hk++hk-I)
,
k-1,2.
A. I: Fedorchenko
152 For the approximation Of gijn
{awldx:,} ijn=D;2)wn,
{dWldXi}ijn=~l(“Wn,
pijn=p
(Wijn)
or
~in=l.i,j=C1((WiZI,I+Wijn)/2),
7
/JY*%,j=
(~i:i,j+~ij")/2.
We similarly approximate d (@w/G’s,) /ax,, d (paw/&) /ax,. With the ratio of adjacent steps 0.8Gz,+‘lh,-‘< 1.3, k=l, 2, this approximation of the second derivatives, being formally of first-order accuracy in hkf 1, will not significantly worsen the accuracy of the solution compared with the case of a regular mesh template (see [lo]). We note that the central approximation used for all the spatial derivatives is especially desirable when investigating acoustic processes in bounded domains without selecting any predominant direction of wave propagation. It can be shown (by the example fk = Ukw, Uk = const), that neglecting g at stage 1) introduces an additional error 0 (zpukdg/as,) , which is fairly small if /3< I, r < 1, into the over-all approximation (4.1). An investigation of the difference scheme in the conditions of the linearized Navier-Stokes equations has shown that on a weakly irregular template this scheme has accuracy 0 ( rm+hZ+$ (I?+) ) (1z=mas (h.12=‘), m = 2 for E = 0.5, M = 1 for e > 0.5 and possesses the necessary local stability criterion r
(rl, r?) , ( l/Re’,
t,--hJ(c+~u~)
(2Ej1.\,
wh,,=l3$‘,
~20.5.
(4.4)
y/Re’Pr’) , c is the speed of sound, and h, = min (h 1?I . h2 k 1).
The quantities Re’, Pr’, c are estimated taking into account the local values of p, T, p, In the numerical solution of a number of essentially non-linear problems the sufficient conditions of stability were close to (4.4). In many real problems for Re > 10 to 102 the value of the time step is mainly constrained by the condition r G 71 (that is, the local numbers K =r (c+ 1u I) /h,,=O (1) ) . which is a necessary physical criterion of solvability of unsteady gas-dynamic processes, where it is important to take into account the mechanism of the propagation of acoustic perturbations. This criterion is practically independent of the form of the difference scheme, explicit or implicit. This scheme generalizes naturally to the case of three spatial coordinates. 2. We will discuss in more detail some features of finite difference methods of investigating acoustic processes, As a rule the dissipative and dispersive properties of difference schemes of gas-dynamics are analyzed by the example of the Cauchy problem for a fust-order hyperbolic system of equations, which is a simplified analog of Euler’s equations. The propagation of plane waves in a quiescent medium is usually considered.
Unsteady subsonic frows of
viscousgas
153
For example, a linear analysis of the difference scheme (4.2), (4.3) shows that it possesses scheme viscosity which exerts a damping action on the longitudinal (along the wave propagation vector) distributions of the parameters, that is, it is analogous to volume absorption. By considering in conditions of the the propagation of a plane harmonic wave &‘-wO esp [i (k,s-w,t)] difference approximation of the linear transport equation &cl&=--dwldx, it is possible for a specified wavelength An to estimate (assuming that h,lh,< 1, T/h,< 1) as in the approach in [ 1 I]) the mean dimensionless coefficient of volume expansion (a,) h (w -esp [- (a,) LX]xesp [ - ( CY,)ht] ) due to the scheme viscosity:
It is obvious that increase of the parameter E > 0.5 enhances the stabilizing properties of the scheme (the damping of high-frequency components of the spectrum), which is sometimes desirable when computing normal discontinuities, flows with high Re numbers etc. The average relative decrease of the phase velocities (that is, the effect of scheme dispersion) is determined by the relation (c,)~z~-- (2rrh,lh,)“, ~~[0.5,11. However such investigations do not take into account the principal features which arise in the numerical solution of the boundary value problems of the dynamics of a viscous gas. Allowances for such factors as the dissipation and dispersion along walls, the acoustic properties of the boundaries, the existence of inhomogeneities of the flow etc. can radically change the approach to the estimates of the scheme errors. Let us assume that in the bounded domain G the spatial modes (independent within the limits of the linearized problem) of the acoustic oscillations are investigated, each mode numbered m (m = 1,2,3 according to the number of coordinates) forming a spectrum of harmonic oscillations with frequencies o,, and wavelengths A,,, n = 1, 2, . . . . Then for each cm, n> component of the spectrum the necessary solvability criteria are connected with the requirement of definite orders of smallness of the quantities di:
(4.5)
Here hy, hw are the characteristic distances between the mesh nodes inside the domain G (along the direction of wave propagation) and close to the boundary I’ (along the normal); 6 011, 6(s) are the thicknesses of the (viscous and temperature) acoustic boundary layers, see [6] ; (c,,)~ is the local phase velocity of wave propagation within the limits of the initial system of differential equations; (c,,)h is the phase velocity in the conditions of the finite-difference equations (that is, allowing for scheme dispersion). The effective attenuation (acceleration) factor (cx,,)~ (the amplitude fmn - exp (-omnt)) can be provisionally represented as the sum (a,,)
p
(a,+a,+arJmn,
(4.6)
154
A. T. Fedorchenko
where 0~~takes into account the losses (sources) in the volume G, (Y,,,allows for the losses close to r due to viscosity and heat conduction, CQ allows for the total flows of acoustic energy through the permeable segments of the boundary (for short ducts, situations are realistic where ICY~ I% ICY”t aw I). The quantity (arnn)h defines the total contribution (in the entire domain ch) of the scheme effects in the attenuation process For example, for the propagation of longitudinal (plane) waves in a circular tube of constant cross section with impermeable walls (in the absence of a mean flow and with weak attenuation [7,
121) 6n(L)/r0~6~) lr,=0(2/o,
Re,)‘“
(c,),/~,~l-6,(~‘/2r~Cl,
(a,) .-0
(o./Re,)
(a,) .=O
Pr=1,
con=6.hnro/cu,
(4.7)
(o,*AW ,
“‘,
where usually in the low-frequency part of the spectrum (a,), of volume absorption due to /J # 0, x+0).
% (cu”)~(here (cx,,)~is the coefficient
In general cases (taking into account the vortex flow, non-linear effects etc.) checking that (4.5) is satisfied may be a separate extremely complex problem.
5. Examples of numerical solutions 1. Using the complete system of equations describing the axisymmetric flows of a viscous perfect gas (the form of the initial system is similar to [3]), a number of problems were solved in the rectangular domain G: {ZX [0, L] , r= [0, I], L=8}. In all the problems the following was assumed: cc= x-i, 7 = 1.4, Pr = 1: the initial field was specified in the form u (x, r, 0) = 0, ~(x,r,O)=p,=l/y,T(x,r,O)=l,x,rE~;onthesidewallu(x, l,r)=O,T(x, l,t)=l, r E .Q = (0, fk), x E [0, L]; on the axis (r = 0, x E [0, L]) we have the conditions of symmetry (1.2). All the calculations were carried out up to fk = 100 to 150. In the domain G meshes were constructed with practically uniform step hX (in the different versions h, = 0.2 to 0.05) and essentially non-uniform step h, (along the wall (/~~)~b * 0.05 to 0.01). The total number of mesh nodes in cj, was in the range 600 - 3200. A repeatedly specified retuning of the mesh could be carried out during the solution by a special algorithm. The calculations (on the BESM-6 computer) were carried out on the basis of the scheme (4.2), (4.3) with v 6 1, E = 0.5 to 1. The time step was specified (from the condition 7 I; 71 in (4.4)) within the limits 0.01 - 0.03. Estimates in accordance with Section 4 subsection 2 showed that a fairly wide band of the spectrum of the longitudinal oscillations wsll from n L 10 to 20) and a narrower band of radial oscillations (w,, from n 5 4 to 8) were resolved without substantial errors in the calculations carried out. 2. In the first series of calculations (Re, = 102, 103, 104) in the statement of the boundary conditions on rl, r2 the model {m, p} = {,%, (0)) S,(L) } was used: pu, (0, r, t) =0.2 (l-r’) T(0,
r, t) =I,
[ l-exp
(-0.2t)]
20,
~~(0, r, t) =O;
(5.1) p
CL, r, t>=p,=lly;
T(L,
r, t)=1,
uz(L,
r, t) =0
if u,CO on r2,
rE [0, I].
155
Unsteady subsonic flows of viscousgas
As a result of the specification on I’1 of a monotonic growth of the flow rate the general development of the flow in the duct was accompanied by the generation of intensive longitudinal anharmonic oscillations with a clearly defined dominant frequency ox. 1 = 0.186, close to the fun&mental characteristic frequency of the one-dimensional quarter-wave resonator w,lo = n/2L = 0.196. Also resolved in the calculations were the radial oscillations (that is, normal waves of higher for orders [ 1,6]) of considerably less intensity ( max 1p (x, 1, t) -p (5, 0, t) 1/p,G10s3 r Z IO) with fundamental frequency w,.r = 3.6 % wxl. We note that the generation of radial oscillations arises not only from inhomogeneous conditions on rl (if dpu, (0, r, t) /drSO, then dp (0, r, t) /6’r+O), but also due to refraction effects in the vortex flow [ 11. For Re, = 102 the attenuation of the oscillations and the establishment of a steady state of the flow with p (0, r, r)/p, = 1.09 was observed for r 2 100. In all the cross sections the steady state profiles of ~1, pul were close to parabolas.
FIG. 1 ,
FIG. 2
Figure 1 shows graphs of m, (t ) /mar m2(t)lmo (m,=m,(~)=m~(m) =x./10). igs l-3 respectively) for Re, = 103, t < 60. We note that the curves of ~(0, 0, G/pa P p (0,0, t)/pa in Figs. 1-3 are averaged over small radial oscillations. For Re, = 102 in the given problem ox,,= (zlL.)=,,, (a,) X,,B (a,) Iti, at least for wxn with n 6 10. As Re, -+ 00(II + 0) for all frequencies wxn we must have low + oVI + 0, Q + or (see (4.6) (4.7)). Considering in G 2 homogeneous flow of an ideal gas with mean velocity (U I )a = MU z 0.1 and putting Qr = --au, Q2 = =(see (3.1) (3.5)) it is easy to show that or = 0 (@,/L) = 0 (10-z). Here the value of or, to a first approximation independent of frequency, is determined mainly by the convective transfer of acoustic energy through the cross section r2, From the numerical solution of the problem for the fundamental tone wX1 it is possible to determine the effective attenuation factors (YZ 0.6, a = 0.03, (Ya 0.02 (for Re, = 102, 103, 104 respectively), which agree with the estimates made. For Re, = 104 (Re, = RecMmax = 2000) an unsteady separation of the flow from the side wall, the development of secondary processes on instability in the viscous shear layer being torn off etc. were observed. These precritical phenomena are discussed individually in [2,3].
156
A. T. Fedorchenko
3. Also solved was a version of the problem (for Re, = 103) in which only one change in (5.1) was made: instead of specifying on I’1 the temperature distribution T(0, r, t) = 1 a condition of the type (2.6) is introduced: T(0,
whereg(r)=
r, t)=l+g(r)
1 forrE
[ (~(0,
[0, l-S],g(r)=
r. t)lp,)(i-i)‘Y-l],
1 - [(r-
rE[O,
1 t&)/S]2 forrE
11,
t4,
[l -6, 11.
This condition is equivalent to specifying a quasi-stationary entropy profile on the segment rl: (x = 0,r E [0,1 -61). In the calculations 6 = 0.2 was used (that is, 6=0 (6,“‘) =0 (a,‘“‘), see (4.7)). As is to be expected the amplitudes of p (0, 0, t)/~~ exceeded only slightly (by not more than 10%) the corresponding values of the version (5.1), the general picture of the course of the process being practically similar. 4. In another version, Re, = 103, in relation to (5.1) only the condition on r2 was changed in accordance with the model (3.4): P c-k I’,
where@(t)=OfortE
t) =Pl+p (t) [m(t) --m* (t) 1,
rE[O,
I],
t4,
[0,tb],P(t)~00.28fort>tb.
The instant ~b = 11.7 corresponds to the first equalization of the outflows ml and m2 in Fig. 1. Therefore, by (3.9, (3.6), on r2 in the linear approximation we have s22 = 00(92 = -1) for t < &, Q2 = 1 (172 14 1) for t > tb, that is, for c > cb there is practically no reflection of the quasiplane wave from the cross section r2, which must lead to the rapid attenuation of the longitudinal oscillations in G. Figure 2 for t Q 60 shows graphs of m 1 (t)/ma, m2 (t)/ma, p (0, 0, t)/pa, p (15, 0, t)/po (curves l-4 respectively). It is seen that for the solution of the non-linear problem the results of the linear estimates for Ql, Q2 are confirmed. The attenuation of the oscillations and the establishment of the steady state of the flow is observed in practice for f 2 50. It is obvious that such an effective attenuation of the longitudkl oscillations (in the entire frequency spectrum computed) could not have been accomplished by introducing some model of volume absorption. 5. A similar attenuating effect was obtained on solving the problem with Re, = 103 also, where unlike (5.1), on rl a model S,, of the form (2.7), (2.8) was used: Pf
(0, r, t>=cp(r) be”(t) --l-f CL)1 t
ph(O,
r, t) =O,
T(0,
r, t)=l
if
u1 (0, r, t) 20,
(5.2) peo(t)=pa{l+0.2[1-exp
(-0.2t)]},
cp(r)=0.2(1-rZ)/0.191p,~=y(1-r’),
p”(t)=p(O,
0, t),
rE[O,
t=Q.
I],
From (5.2). (3.3) we can estimate Qr = 0 (-I), In1 I < 1.~2 * - 1. that is. in this case also rapid attenuation of all the frequencies w,, can be expected. Figure 3 shows curves 1-4 corresponding to m, (t) lm,, m2 (t) lm,, p (0, 0, t) lp,,, Establishment of the steady state (with the outflow m, = ml = m2 differing p,“(t) /pa. negligibly for the parameters of (5.2) from m, in problem (5.1)) is observed for t>,50.
1.57
Unsteady subsonic flows of viscous gas
'a
20
4D
t
&!7
FIG. 3
6. Using the model {p, I)} the problem of the acceleration of a subsonic flow of an originally quiescent gas when there is a monotonic increase of the pressure on rl was solved numerically: ~(0, r. t)=p,{l+d[l-esp
(--h-t)]}.
T(0,
r, t)=l,
~~(0. r, t) =O, p(L,
ford>O,k>O;ul
r, t) =pa=const
>,Oonrl,r2.
We recall that the solution of the linear problem of the acceleration of the flow of a viscous incompressible fluid in an infinitely long circular tube (on the instantaneous application of a specified longitudinal pressure gradient to the quiescent fluid) was first given in [ 131, and later (and apparently independently) in [ 141. Under real conditions the accelerated flow of a compressible gas in a duct of finite length is usually realized by varying the pressure at one end of the duct, which may cause the development of longitudinal acoustic oscillations of considerable amplitude. Here the solution of the problem for Re, = 100, d = 0.1, k = 0.2 is shown. The choice of such a low Re, number was determined by the fact that for these values of Re,, L the effects of viscosity show themselves clearly even in the course of a comparatively short time interval. The calculation of this version was continued up to ff = Re, = 100, when in practice the attenuation of the oscillations and the establishment of the steady state of the flow was observed. If we introduce the relative time t* = t/Re,, independent of the choice of the characteristic velocity co, then the total “acceleration time” tf* x 1 agrees with the solution of [ 141.
Figure 4 shows graphs of w,(t), u:,(i) (w,(t) =mj(t)/n), ~(0, T’, t)/pa, p(L/& 0, t)/p, (curves l-4 respectively). For comparison curve 5 gives the graph of wf(r*)wl (lOO)/wf (I), where wf (t*) is the outflow from the solution of [ 141. If we estimate the “switching in” time of the constant longitudinal pressure gradient of value tR * I/k * 5. then by a shift of curve 5 to the right by At * 5 it will be close to the variation of the “average” outflow (WI + ~2 )/I!.
From the curves of Fig. 4 it is possible to determine approximately for the fundamental frequency ox 1 * n/L the characteristic attenuation time tl =Z15 to 20, and also the decrease of the phase velocity (cl = 0.9), which agrees with the estimates by formulas (4.6) (4.7) (cL+O, (a,) ix l/t,=1/16> (a,) ,).
158
A. T. Fedorchenko
1
P/Pa
FIG. 4
FIG. 5 Figure 5 shows the profiles of the longitudinal velocity ~1 (0, r, t), u1 (L, r, t) at various instants of time. It is obvious that the effect of friction at the wall, as in the solutions of [ 141, is gradually transmitted along the radius to the axis. For t = 100 the profiles of ~1, pul in all cross sections are almost parabolic in shape, WI = w2 = 0.112, I (WI - wz)/wl I S 0.3%. We note that in this problem in conditions of quasi-parallel flow ( I~U~lic~]]~Zllo, (IdplLk& >~jQddyIIG) the formation of the profiles UI (r) (including that on I’l) is primarily determined by the terms apldx, (r Re,) -9 ( prau,ldr) ldr in the equation of motion along the Ox axis. It is characteristic that the asymptotic solutions (as t + -) of all the problems of subsections 2-6 define on a finite segment of the cylindrical duct some analogs of the Poiseuille flow for the case of a compressible medium. Therefore, using the models presented we can investigate the nonlinear effects of the action of acoustic oscillations (for example, in the spectrum of the characteristic resonance frequencies w xn, car,, of the domain G) on this important type of flow. The author thanks A. A. Samarskii, K. I. Artamonov and 0. S. Ryzhov for discussing the results and for valuable comments TramWed
by
J. Berry
159
Unsteady subsonic flows of viscousgas REFERENCES 1. NAYFEH, A. H., KAISER, J. E. and TELIONIS, D. P. Acoustics of aircraft engine-duct systems AIAA J., 13,2, 130-153,1975. 2. FEDORCHENKO, A. T. On non-linear acoustic effects in the flow of a viscous gas in a cylindrical duct of finite length. Do&LA&ad Nauk SSSR, 237,6, 1315-1318, 1977. 3. FEDORCHENKO, A. T. On calculations of two-dimensional unsteady flows of a viscous gas in a short cylindrical dust with an end blast. Izrl. Akad. Nauk SSSR. Mekhan. zhtdkosti igaza, No. 1,9-17,
1919.
4. THOMAS, P. D. Boundary conditions for implicit solutions to the compressible Navier-Stokes equations in finite computational domains. AM4 Comput. Fluid Dynamics Conj: Williamsburg Vu, 1979, pp. 14-26. 5. ROZHDENSTVENSKII, B. L. and YANENKO, N. N. Systems of quasilinearequations and their applications to gas dynamics (Sistemy kvazilineinykh uravnenli i ikh prilozheniya k gazovoi dinamike), Nauka, Moscow, 1978. 6. ISAKOVICH, M. A. Generalacoustics (Obshchaya akustika), Nauka, Moscow, 1973. 7. LANDAU, L. D. and LIFSHITS, E. M. Mechanics of continuous media (Mekhanika sploshnykh sed), Gostekhizdat, Moscow, 1953. 8. FEDORCHENKO, A. T. Numerical investigation of some MHD flows of a viscous compressible gas in a rectangular cavity. Izv. A&ad. Nauk SSSR. Mekhan. zhidkosti igaza, No. 5,27-33, 1975. 9. FEDORCHENKO, A. T. Numerical investigation of an acoustic resonance in a rectangular cavity in a viscous gas flow. Akustich zh., 24,5, 75 l-759, 1978. 10. SAMARSKII, A. A. and POPOV, Yu P. The difference schemes of gas dynamics (Raznostnye dinamikil. Nauka, Moscow, 1975. 11. ROACHE, P. J. Computationalfluid dynamics Hermosa PubL, Albuquerque,
skhemy gazovoi
1976.
12. RAYLEIGH. Theory of sound (Teoriya zvuka), VoL 2, Izd-vo in lit, Moscow, 1955. 13. GROMEKA, I. S. On the theory of motion of a f7uid in narrow cylindrical tubes (K teorii dvizheniya zhidkosti v uzkikh tsilindricheskikh trubkakh). Kazan, 1882. 14. SZYMANSKY. P. Quelques solutions exactes des equations de I’hydrodynamique de fluide visqueux dans le cas d’un tube cylindrique. J. Math. pures et appl, ser. 9, 11, 1, 67-107, 1932.
0041-5553/81/050143-17$07.50/O 0 1982. Pergamon Press Ltd.
A METHOD FOR THE NUMERICAL INVESTIGATION OF UNSTEADY SUBSONIC FLOWS OF VISCOUS GAS IN DUCTS* A. T. FEDORCHENKO Moscow (Received 14 February 1980)
A NUMBER of computational models for solving non-linear problems of unsteady subsonic flows of a viscous gas in a finite segment of a duct are presented. A difference scheme is presented for the integration of the complete Navier-Stokes system of equations for a compressible medium. Some general features of the application of mesh methods for investigating three-dimensional acoustic processes in ducts are considered. Examples are given of the realization of computational models and of the difference scheme in a series of unsteady problems of the development of the axisymmetric flow of a viscous gas in a short cylindrical duct.
Introduction At the present time in the general development of computational methods of gas dynamics the method of computing essentially unsteady subsonic flows of a viscous gas in ducts using the multidimensional Navier-Stokes equations is still insufficiently developed. The absence of solutions of non-linear boundary value problems of this class is objectively explained both by the difficulties in constructing models of unsteady subsonic flows of a viscous gas in bounded regions (this is due to specific requirements imposed on the boundary conditions), and also by problems of the correct resolution of the three-dimensional acoustic processes by a difference scheme. In particular, the investigation of such topics would facilitate the solution of important applied problems in non-linear acoustics [ I 1. In this paper a group of models is presented for computing the unsteady subsonic flows of a viscous heat-conducting gas (using the complete three-dimensional Navier-Stokes system of equations) in a finite segment of an arbitrary duct. The acoustic characteristics of penetrable end barriers of various types are investigated, which enables methods of efficient action on the longitudinal oscillations in the duct to be indicated (for example, for the purpose of damping the low-frequeny components of the spectrum). For application to a wide class of gas-dynamic problems a difference scheme for integrating the complete Navier-Stokes system of equations for a compressible medium is proposed. The fundamental characteristics of the scheme are given. Some necessary mesh criteria for the solvability of the acoustic processes in the ducts are determined. As examples of the use of the proposed models and the difference scheme numerical solutions are presented of a series of unsteady axisymmetric problems in which non-linear oscillatory processes in the development of the flow of an originally quiescent viscous gas in a short cylindrical duct are considered. *Zh. vjkhisl. Mat. mat. Fiz., 21,5, 1215-1232,
1981.
143
A. T. Fedorchenko
144
1. General formulation
of the problem
Consider a bounded domain G, which is a finite segment of a duct with rigid impermeable side walls. For the general formulation of the problem we introduce some curvilinear orthogonal coordinate system {OX. Oy, 02). assuming that the axis Ox is directed along a slightly curved nominal axis of the duct. We assume for simplicity that the permeable end boundaries rj, j = 1,2, are elements of the coordinate planes x = xl, x =x2, x2 - x1 = L. The duct wall, the boundary r,- (I’==~,uI’,-ur2, I’,,,=r,, n I-,,, C=GL r) , is defined by some continuous piecewise-smooth surface _!Z(2, zj, z) =O. The subsonic unsteady flows of a viscous heat-conducting gas are modelled numerically in the domain G. Without significant loss of generality we wilJ use the model of a one-component ideal gas, governed by the equation of state p = pT/-y, y = cPo/c,o = const, which is in equilibrium in the conditions of the low-frequency approximation. Then, for the five independent parameters of the flow ~1, ~2, ~3, p, T the initial system of five differential equations includes the continuity equation, three equations of motion (the complete Navier-Stokes equations) and the energy equation. These equations are complemented by the local relations p = p (T, p), X=X (T, p), . . . . Here and below we consider the dimensionless values of the coordinates (x, y, z), the velocities (~1, ~42,ug), the time r, the density p, the pressure p, the temperature T, the coefficients of shear viscosity p (the coefficient of volume viscosity is ignored) and the thermal conductivity x. These quantities are expressed in the units rg, co, r&g, po, poco2, To, ~0, x0 , respectively. The characteristic values chosen were rg the characteristic half-width (radius) of the duct, To the characteristic temperature of the walls, co (To) the adiabatic velocity of sound, po, ~0, x0, cP,, , the density, viscosity, thermal conductivity and specific heat at the temperature To and some characteristic (initial) pressure in the duct. We define the fundamental similarity criteria:
In the domain G for the interval R = (0, tk) we can formulate the unsteady boundary value (more accurately, mixed) problem with the initial field {u, p, 2’)~ at t = 0 and boundary conditions on r. On the impermeable wall of the duct rw, conditions of the following form are usual:
u=o,
a,dTldn-I-a?T=a,, ai=ai(r, 1).
d,,
t4.
(1-l)
Some element of the boundary may be a plane (axis) of symmetry, where the conditions
u,,=o,
du,ldn=dpldn=dTldn=O,
(1.2)
are specified, and u,, U, are respectively the normal and tangential components of the velocity vector. The specific models of the flow in the domain G discussed below are primarily defined by the method of specifying the boundary conditions on rl, r2.
Unsteady subsonic flows of viscousgas
145
2. Basic computational models The proposed non-linear models of unsteady subsonic flows in the domain G include as important elements two types of permeable boundary surfaces (we denote them by S,,,*, S,*), intended for specifying the conditions on ri. To be specific we assume that at all the points r E T’i considered we have U, = (u, n) # 0. The upper suffvres “t”, “-” denote the intake of the flow at the given point T’i: (u, n) > 0, and the outflow: (u, n) < 0, n is the inward normal. 1. Let us first consider the &-surface (some simplified versions of this model were proposed previously in [2,3]). By the general model of S,, at each point r E S, + C Pi, t E 52, conditions are used of the form
\‘(r. f)=f(r.
t, p(r. t)).
7’!r, t) =g(r.
f-{:i+
(2.1)
!:> in}.
t, p(r, t)),
(2.2)
where f, g are given functions. Here and below Vi = pF ui (6 is a specified real parameter, in general different for i = 1, 2, 3). Sometimes instead of (2.1) it is preferable to specify the function $ and the unit vector b: P’luI =$(r,
L p(r, t)),
u=lulb(r,
t).
(2.3)
In the case U, < 0 (on S,,, - ) a single condition is imposed of the form
t’Tl(r, t)=(V.
n)=f),(r,
t, p(r,
t)),
rES,, -=I?,,
J&.
(2.4)
We will discuss separately conditions of the type (2.2). In the simplest case, ag/ap = 0, the temperature distribution
I’r=g (r, t) ,
rEr,+,
t4.
(2.5)
can be specified on ri+. Physically this condition corresponds not to the cross section of the developed flow, but rather to a permeable wall (lattice), and therefore it would be more suitable for a permeable segment of the boundary T’,+ Obviously, even for the whole boundary Pi = Pi+ it is possible to use the more general relation T(r, t)=gl(r,
L) [P(r,t)]‘:-‘“‘+g?(r,
t),
glao,
(2.6)
which may be equivalent to the specification of a slowly varying entropy profile s (r, t) over a large part of the cross section Pi+ (except for a neighbourhood close to Pwi with characteristic width of the order of the thickness of the acoustic boundary layer 6, = max (SF’ ,a%’ > along the wall of the duct, see (4.7)). We note that in a number of cases the use of (2.5) instead of (2.6) should not lead to significant differences in the solutions, at least for the low-frequency band of oscillations (see below remark 2 and subsection 3 in section .5)., We present some particular cases of the model S,.
146
A. T. Fedorchenko
Within the limits of the S ,+,-model the instantaneous values of V, (r, t) at the point r E F’j are determined (with possible time delays ri) by the drop between the pressure p (on the “mternal” side of I’i) and the “external”
pressure pe:
v,, (r, t) = cp(r, t) I AP Ik sign (APL AP =
:i
s:,
me (r, tl) -
f2 (r, t, VI (r, t)),
vt (r, t) =
T (r, t) = g (r, t,
i
where r,>O,
920,
p,>O,
0, t),
i=l,
2,
V3 = fs (r, 4 VI (r, th
(2.7)
P Cry W;
t) 1bplk sign (Ap),
V,,(r, t) =q(r,
Sm,-:
k > i!i=t-Ti(r,
P (1, tz),
g>O, and
In some cases the distributions
f2, f3 are specified functions.
p (r, t), pe (r, t) can be replaced by the average values po (t),
p,o (t) over the cross section ri:
v, (r, t) =(r
(r, 0 I APO Iksign
(Ap”),
ApO=p.O(t,)
... .
-p”(t,),
(2.8)
In the simplest version of the model S, (we denote it provisionally by S,,, 0) it is assumed that the vector V at the given point r E rl .+ is specified, that is, it is independent of the oscillations of the pressure inside the domain G: Smo+: V (r, I) =f (r, t),
T(r,
t),
Smo-: V,(r, t) =fn(r,
t)=gk,
6 p(r, t)), (2.9)
tEQ2.
rEl?j,
The model Snto+ can be formally regarded as the result of the passage to the limit in Sm+,+: pe + 00, cp+ 0, p < const. However for S, 0 - such an interpretation
is invalid (pe 2 p for (u, n) < 0).
2. By the model S,, for r E Pi, r E R there may be conditions
T(r,t)=g(r,t,p(r,t)),
~(r,t)=f,(r,t,V,(r,t)),
S,+. -1
&-:
of the form
u/lul=b(r,t): p (r, t)
(2.10)
=fP b, 4 Vn(r, t) ) ;
fp, b, g are specified functions. In a particular case (the model S,n) af,/aV, = 0. In &other
particular
S,+: p(r, SP-:
p(r,
version of the Sp-model the conditions
t)=fr(t, l) =fP(t,
mj(tj)), W(tj)
T(r, )7
t)=g(r,
r=rf,
4 p(r,
t)),
b=b(r,t),
(2.11)
t42,
are possible, where fi = t - rj (t), ~j > 0 is a specified delay parameter,
rni is the instantaneous mass outflow in the cross section I’i (in calculations for ri = 0 the approximate values of rnjdetermined in the cross sections x1 + hX or x2 - hX, 0 < hx < 1 may be used).
In some cases (for example, if a~p+azV’,=a,, equivalent
to S,.
ai=ar(r,
t) , rd?,)
the model SP may be
Unsteady subsonic flows of viscousgas
147
The models S,, S, can be naturally generalized and made more complicated by introducing “controlling boundary conditions” of the type
V(r, t)=f(r,
t
)...,
Jr),
rdm+crj,
Ed&
(2.12)
where W’ is a finite set of parameters of the process investigated from the solution of the problem on the point set G’XQ’, G’cG, R’c [t,, t) d2. Similar analogs are possible for fP, fn,g,b, . .. also. In particular, conditions (2.7), (2.11) with ri > 0 can be reduced to the simple form (2.12). In general, the process of numerical modelling of the gas-dynamic phenomena in the domain C inescapably involves the need for a continuous investigation in the computational process of some set of parameters of the process of solving the problem in the preceding time interval (with local estimates of the number M and the sign of (u, n) in I-j, with an analysis of the spectrum of the acoustic oscillations generated and selection of the next algorithm of operation on the oscillatory process etc.). Remark 1. The assumption of the subsonic nature of the flow in c (in every case, M < 1 close to I’) imposes certain constraints on the functions f, fp, g, . . . in (2.1) - (2.11). In some simple cases (see section 5) preliminary estimates of the range of parameters in the solution of the problem are possible. If it is impossible to provide such preliminary estimates the tracking of M and also of other parameters in G may be accomplished by using controlling conditions of the type (2.12). 4.
An approximate analysis of the initial system of differential equations close to rj enables the number of independent boundary conditions in the models S, +, SP* to be validated. Indeed, let the boundary Fj be a normal cross section of a quasi-parallel subsonic flow. Then near each point rj the initial system of five equations can be interpreted as locally hyperbolic in the variables x, t (see, for example, [4]). Here, all the gradients along the axis Ox are assumed to be small, and the terms with derivatives along Oy, Oz are interpreted as distributed sources or force functions. A subsequent characteristic analysis [4, 51 of the system leads to the result that the required number of boundary conditions at each point r E I’j+ must be N+ = 4, and at the point r E rj- it will be N- = 1 (correspondingly, for the two-dimensional problem in x, y coordinates we have N+ = 3, N- = 1). It must be stipulated that a similar result is invalid for estimating the number of boundary conditions on a rigid permeable wall of the duct rw. In this case with the dominating role of the parabolic processes for ~2, ~43, T (with respect to x, f in the equations of motion and energy) the number of boundary conditions may amount to N+ = N- = 4, although here also the basic elements of the models S,, SP can be used to determine the connection between ul, p, T. It is important that the very nature of the proposed system of boundary conditions (2.1) (2.11) differs significantly from the traditional approaches of various authors (for example, in [4] for I?,+ it is proposed to specify the angles of approach of the flow, the stagnation temperature and the stagnation pressure, that is, conditions obviously inadmissible for the modelling of viscous subsonic flows, especially for taking into account non-linear acoustic processes).
148
A. T. Fedorchenko
5. It is natural that in (2.1) - (2.11) the number of independent boundary conditions is always less than five. Therefore, in the numerical solution of the problem the parameters at the boundary points must be calculated by introducing additional difference algorithms, containing to some degree elements of spatial extrapolation from the internal points of the domain. These algorithms are based, in particular, on the characteristic dependence of the pressure at the points r E S, C r. I E f’,+, on the process in the entire domain G by means of the acoustic connection ? (M < 1 in G). For example, when the difference scheme of Section 4 is used the pressure ,;+r (r=S,, r=r,, tn+‘= (n+l) ‘c) can be computed by normal extrapolation of second-order accuracy from internal points of the domain G, already computed on the layer tn+l (sometimes this procedure is carried out by using a difference analog of the continuity equation, written in conservative form for the boundary cells of the mesh). Then conditions (l.l), (2.1), (2.2) enable us to compute VF+l and T{” on S,+, &?+l, $?+I,. . . . Various extrapolation procedures (for example of the type (d2~!dn2)ri+i=O)can be used to compute Tf!‘l on S,,- and S,-. V;+l on S,+. 6. We will determine the general model of subsonic flow in the domain G as a combination of the S, and SP surfaces used to specify the boundary conditions on l?j. Therefore, three basic types of computing models can be introduced: {m, p}={S&& {m,
S&)1,
(2.13)
~}={~n(~i)~Sm(~A))~
{P,P}={S,(5i),Sp(5k)},
ifk,
i=l,2,
F;=l,2.
We mention briefly some of their features and limitations. If in terms of the model (2.13) the surfaces rj are identified with the end cross sections of a complete length of duct, then in this case it is difficult to take into account local end effects close to I’,+,i(the vortex flow past angles and sharp edges, diffraction wave phenomena etc.), depending on the shape of the boundaries and the conditions of the flow outside the domain G. However, this problem can be solved if, using the same models of S,, S,, we consider the system of “ a duct inside a duct” (that is, a comparatively narrow duct connects domains with a sharp expansion of the transverse cross section). Moreover, the models of (2.13) are correct only for wave processes with weak discontinuities, since the conditions of reflection of shock waves from the boundaries I’i are more complex in nature, going beyond the limits of the models S,,,, S,. Using the model SPo it is practically impossible to specify in the “input” cross section rj (in the sense of the direction of the averaged mass flow) arbitrary profdes Vl with characteristic values of the transverse gradients (especially with a definite value of ml (r)). Therefore, in the majority of cases, except for a number of specific problems (see subsection 6 in Section 5, the model $0 can be usefully used only in the “output” cross section.
3. The acoustic characteristics of the permeable boundaries By linearizing the boundary relations (2.1) - (2.11) it is possible to estimate the acoustic properties of the permeable boundaries r’.
Unsteadysubsonic flows of viscousgas
149
1. Let us consider a small element dTj=rj (to be specific, j = l), on which, using the model (2.9), we put V,=pFui=f (t) , where 6 = const, f(t) is a function, which is fairly smooth and monotonic for t E (r 1, t2)c Cl, and is practically unchanged within a characteristic period of oscillation fp = 2rr/w < t2 - rl . Close to dl?l* (on clTl+ for g2 = 0, gl G=const > 0 in (2.6)) we will consider small normal adiabatic oscillations (u,=u=ii+u’, p=pi-p’, p =pSp’, p’=p’/~‘, / u’iii 1
&i/C=---jii?
sign (li), (3.1)
‘I= (l+q)/(l-q)
for q*1.
It is characteristic that 9, n are real and are independent of the frequency of oscihation. We note that 9 ~-l(~~O)f~r~~~/ii(thatis,1~1>1i;9-+O(~~1)asC;~O.For~E(O,~~/~~)we have 9 < 0,~ < 1 on d rl+ and 9 > 0, n > 1 on d I’1 -. The case j = 2 can be discussed similarly. However, significantly different results are possible when the oscillating functionsf (t) are specified on&O C q.
2. Let us consider the normal harmonic oscillations (U‘=uo’exp (--id), p’==po’ exp (--Lot) ) close to the element d I’l, on which the model (2.7) is used with k = 1, ~1 = 0. We assume that cp(t), pe (t), ~2 (t) are practically unchanged in the course of the characteristic period of oscillation rp. Then for drl*
Putting 11~ il5.
T:=& p’Z”(c. (,LI~-~). we can rewrite (3.2) in the form
E---
q
u 2L.t c [ Pe-P
I
Fz-
Wra+j(pe-_P)) Ipp-jd
(3.3)
*
For example, 9 * 0,~ = 1 for zx-rP/(p&,-I_i). For [= 1 (VI = pul) we have Yp+c (PC-p) =p (y-1) +p20, that is, 9 < 0 on drl*. We note that (3.3) reduces to (3.1) only for the case when (u, n) > 0 as (Y+ 0, pe + m, P < const. Remark 2. Just at the points r E lYlwith the condition (2.5) for ag / dr=Owe have p’=:yp’ / T?. The relation p’=p’ / ~2 is true only for rs.r,t 6:’ (see (4.7)), which by the way does not change the order of the estimates (3.1) - (3.3). In any case the portion of the acoustic energy absorbed in such a case due to heat exchange (see [6,7]), can be estimated by the quantity AIU/ W= 0 (0, / R&Pr)‘+ usually extremely small in the low-frequency band.
The total effect of the weakly curved surface J’i (with the distributed quantity 9 (r, t), I E T’i)can be estimated approximately by “averaging over the cross section” the values of the conductance Qj and the reflection coefficient nj:
Qj(t> m
+JJq(r,t)da-(Qj),,+i(Qj)~~, ’ r1
j=1,2,
A. T. Fedorchenko
150
= 0, then practically total transmission of the quasi-plane normal wave through the If CQj>im cross section rj (that is, I qj I < 1) must be observed if (Ql Xe = -1 or (Q2)re = 1. Putting (Ql),, > 0 or (Q2)re < 0 may lead to the amplification of the longitudinal (along Ox) oscillations in the domain C. 3. Let us now consider small normal oscillations close to the boundary r2, on which by (2.11) we have put Pb2,
Y, 2, t> =fp
0) i-P@) i.m (0 -6
t’=t--tm,
(t’> I,
Sm),O,
(3.4)
where
The quantity E2 can sometimes be estimated as follows: if on rl the model Smo+ with ~=l,&n~/at~Oisused,thenE~=:m~,andifonI’~ themodel,Spisused,thaniii2~(ml +mz) 12, and so on. Let us assume that fP, 13are slowly varying specified functions (that is, practically constant within a period), and the gradients aTjay, aT/az are small close to F2. Then if rm = 0 then for the quantities uO’, pO’, pO’ averaged over the cross section r2 (with area 02) it is easy to obtain p’=@?[ u,‘p,+p,‘U,]. Putting pG’~p0’,k2, B=?cT?, we have p,‘c, (I-B&l&*)
=Buo’Po&
(3.5)
Q+ (if,--Bii,/&) lB= [ &-B.U, sign (U,) l/B.
(3.6)
We note thst Q2 = 0 for B~:Co2/ii,, U,+O. For C, = 1, H, < 1 the quantity Q2 = 1 (that is, 1~2I Q 1) is attained for B = 1. By specifying the parameter r m > 0 in (3.4) we can obtain various complex values of Q2 (t, w).
4. The difference scheme 1. For the numerical integration of the initial system of differential equations we can use an explicit difference scheme of the predictor-corrector class with the property of over-all approximation, which is a generalization of the scheme previously proposed in [ 81, This scheme, satisfying a number of specific requirements, was successfully realized in an entire series of essentially unsteady problems of the dynamics of a viscous gas [Z, 3,9]. We note that, in general, the models of Section 2 are independent of the type of numerical scheme used. We briefly present the main features of the difference algorithm by the example of the model equation for the vector function w (x, t), x = (xl, x2) E G, t E 52: (4.1)
fR=f(w),
P.=P (w> f
~=O(l/Re)
k, l=l,
2.
Unsteady subsonic flows of
gas
W’SCOUS
151
We will consider some unsteady boundary value problem in the domain G with boundary I’ consisting of rectilinear segments parallel to the coordinate axes. We will construct a non-uniform rectangular mesh (formed from straight lines parallel to the coordinate axes) and consider the set Ch of nodes (z,‘, x2j) EC, eh=GhU l?,,. We will also introduce the set Gh * = Cl * U G2 * of “complementary” nodes situated at the middles of the segments of the coordinate straight lines connecting the nodes of the basic mesh: Gi’= {iPhJ} = {xi=
(~i’+x:*~
Gz*= {s.+*“‘} = {zi=z:,
) /2,
IC~=Z~~, (xi, s,) EG},
zz= (s:+zf’
) /2,
(CC,,x,) EG} .
We define the time mesh: &G
p+o= (n+a)T;
{&O,
cJ”E, 1; n=O, 1,2,.
We will denote the mesh functions at the nodes ChXor’,
..},
&E[0.5, 11.
GhXure by w, w* respectively.
Then the difference scheme for (4.1) is realized within an irregular 9-point template of ~??hin two consecutive stages: 1) at the predictor stage W;,~~~~j=(Wij”+W~*,,j)/2--EZ[DI:)’,,+(~~2’
(wjn+W:j*i)12
Wr.L=
+Di’,“:
/2+0;3f”,
--ET[ (is,“‘+S,‘::) Dj’:I,,f”=*
DI:)l,,fn=~(f,=,,j-f~jn)lhl*‘,
~~‘)fb(p
+ai+'(ftF+!,j-fij")
(f
~i(“fn=az-‘(fu”-f*;-l)
)/Z]f”,
lh,*‘,
(fC,j+l-fzl”)
(4.2)
7
+az+‘(fi;+i-fij”).
h,*’ = I2,‘~-3i
tt1
hl*L Id-& kl I,
1,
ah*‘=hk=‘/hk” (h,+‘+h,-‘) ;
2) at the corrector stage
n+1
Wij
=wij “-t[
Dd“f=2a,-’
(jT*-f&,,)
Do(2)f=bz-i(fz-fi,j-~/a>
fi=Vfijn+
=[O,
+2a,+i(fi;~/z,j-~~)1 +2~+‘(fi~j+~/z-_T*)
(3-V)
fz=Vfijn+(l-V)
(D~‘)+D0(‘))f-~gij”],
)
(4.3)
(b2-‘f~,j-l,*+b2+‘fi~]+1!I)1
(bi-‘f:-‘h,j+bi+‘fi+lil,j),
II,
t&*l=h=‘/ (hk++hk-I)
,
k-1,2.
A. I: Fedorchenko
152 For the approximation Of gijn
{awldx:,} ijn=D;2)wn,
{dWldXi}ijn=~l(“Wn,
pijn=p
(Wijn)
or
~in=l.i,j=C1((WiZI,I+Wijn)/2),
7
/JY*%,j=
(~i:i,j+~ij")/2.
We similarly approximate d (@w/G’s,) /ax,, d (paw/&) /ax,. With the ratio of adjacent steps 0.8Gz,+‘lh,-‘< 1.3, k=l, 2, this approximation of the second derivatives, being formally of first-order accuracy in hkf 1, will not significantly worsen the accuracy of the solution compared with the case of a regular mesh template (see [lo]). We note that the central approximation used for all the spatial derivatives is especially desirable when investigating acoustic processes in bounded domains without selecting any predominant direction of wave propagation. It can be shown (by the example fk = Ukw, Uk = const), that neglecting g at stage 1) introduces an additional error 0 (zpukdg/as,) , which is fairly small if /3< I, r < 1, into the over-all approximation (4.1). An investigation of the difference scheme in the conditions of the linearized Navier-Stokes equations has shown that on a weakly irregular template this scheme has accuracy 0 ( rm+hZ+$ (I?+) ) (1z=mas (h.12=‘), m = 2 for E = 0.5, M = 1 for e > 0.5 and possesses the necessary local stability criterion r
(rl, r?) , ( l/Re’,
t,--hJ(c+~u~)
(2Ej1.\,
wh,,=l3$‘,
~20.5.
(4.4)
y/Re’Pr’) , c is the speed of sound, and h, = min (h 1?I . h2 k 1).
The quantities Re’, Pr’, c are estimated taking into account the local values of p, T, p, In the numerical solution of a number of essentially non-linear problems the sufficient conditions of stability were close to (4.4). In many real problems for Re > 10 to 102 the value of the time step is mainly constrained by the condition r G 71 (that is, the local numbers K =r (c+ 1u I) /h,,=O (1) ) . which is a necessary physical criterion of solvability of unsteady gas-dynamic processes, where it is important to take into account the mechanism of the propagation of acoustic perturbations. This criterion is practically independent of the form of the difference scheme, explicit or implicit. This scheme generalizes naturally to the case of three spatial coordinates. 2. We will discuss in more detail some features of finite difference methods of investigating acoustic processes, As a rule the dissipative and dispersive properties of difference schemes of gas-dynamics are analyzed by the example of the Cauchy problem for a fust-order hyperbolic system of equations, which is a simplified analog of Euler’s equations. The propagation of plane waves in a quiescent medium is usually considered.
Unsteady subsonic frows of
viscousgas
153
For example, a linear analysis of the difference scheme (4.2), (4.3) shows that it possesses scheme viscosity which exerts a damping action on the longitudinal (along the wave propagation vector) distributions of the parameters, that is, it is analogous to volume absorption. By considering in conditions of the the propagation of a plane harmonic wave &‘-wO esp [i (k,s-w,t)] difference approximation of the linear transport equation &cl&=--dwldx, it is possible for a specified wavelength An to estimate (assuming that h,lh,< 1, T/h,< 1) as in the approach in [ 1 I]) the mean dimensionless coefficient of volume expansion (a,) h (w -esp [- (a,) LX]xesp [ - ( CY,)ht] ) due to the scheme viscosity:
It is obvious that increase of the parameter E > 0.5 enhances the stabilizing properties of the scheme (the damping of high-frequency components of the spectrum), which is sometimes desirable when computing normal discontinuities, flows with high Re numbers etc. The average relative decrease of the phase velocities (that is, the effect of scheme dispersion) is determined by the relation (c,)~z~-- (2rrh,lh,)“, ~~[0.5,11. However such investigations do not take into account the principal features which arise in the numerical solution of the boundary value problems of the dynamics of a viscous gas. Allowances for such factors as the dissipation and dispersion along walls, the acoustic properties of the boundaries, the existence of inhomogeneities of the flow etc. can radically change the approach to the estimates of the scheme errors. Let us assume that in the bounded domain G the spatial modes (independent within the limits of the linearized problem) of the acoustic oscillations are investigated, each mode numbered m (m = 1,2,3 according to the number of coordinates) forming a spectrum of harmonic oscillations with frequencies o,, and wavelengths A,,, n = 1, 2, . . . . Then for each cm, n> component of the spectrum the necessary solvability criteria are connected with the requirement of definite orders of smallness of the quantities di:
(4.5)
Here hy, hw are the characteristic distances between the mesh nodes inside the domain G (along the direction of wave propagation) and close to the boundary I’ (along the normal); 6 011, 6(s) are the thicknesses of the (viscous and temperature) acoustic boundary layers, see [6] ; (c,,)~ is the local phase velocity of wave propagation within the limits of the initial system of differential equations; (c,,)h is the phase velocity in the conditions of the finite-difference equations (that is, allowing for scheme dispersion). The effective attenuation (acceleration) factor (cx,,)~ (the amplitude fmn - exp (-omnt)) can be provisionally represented as the sum (a,,)
p
(a,+a,+arJmn,
(4.6)
154
A. T. Fedorchenko
where 0~~takes into account the losses (sources) in the volume G, (Y,,,allows for the losses close to r due to viscosity and heat conduction, CQ allows for the total flows of acoustic energy through the permeable segments of the boundary (for short ducts, situations are realistic where ICY~ I% ICY”t aw I). The quantity (arnn)h defines the total contribution (in the entire domain ch) of the scheme effects in the attenuation process For example, for the propagation of longitudinal (plane) waves in a circular tube of constant cross section with impermeable walls (in the absence of a mean flow and with weak attenuation [7,
121) 6n(L)/r0~6~) lr,=0(2/o,
Re,)‘“
(c,),/~,~l-6,(~‘/2r~Cl,
(a,) .-0
(o./Re,)
(a,) .=O
Pr=1,
con=6.hnro/cu,
(4.7)
(o,*AW ,
“‘,
where usually in the low-frequency part of the spectrum (a,), of volume absorption due to /J # 0, x+0).
% (cu”)~(here (cx,,)~is the coefficient
In general cases (taking into account the vortex flow, non-linear effects etc.) checking that (4.5) is satisfied may be a separate extremely complex problem.
5. Examples of numerical solutions 1. Using the complete system of equations describing the axisymmetric flows of a viscous perfect gas (the form of the initial system is similar to [3]), a number of problems were solved in the rectangular domain G: {ZX [0, L] , r= [0, I], L=8}. In all the problems the following was assumed: cc= x-i, 7 = 1.4, Pr = 1: the initial field was specified in the form u (x, r, 0) = 0, ~(x,r,O)=p,=l/y,T(x,r,O)=l,x,rE~;onthesidewallu(x, l,r)=O,T(x, l,t)=l, r E .Q = (0, fk), x E [0, L]; on the axis (r = 0, x E [0, L]) we have the conditions of symmetry (1.2). All the calculations were carried out up to fk = 100 to 150. In the domain G meshes were constructed with practically uniform step hX (in the different versions h, = 0.2 to 0.05) and essentially non-uniform step h, (along the wall (/~~)~b * 0.05 to 0.01). The total number of mesh nodes in cj, was in the range 600 - 3200. A repeatedly specified retuning of the mesh could be carried out during the solution by a special algorithm. The calculations (on the BESM-6 computer) were carried out on the basis of the scheme (4.2), (4.3) with v 6 1, E = 0.5 to 1. The time step was specified (from the condition 7 I; 71 in (4.4)) within the limits 0.01 - 0.03. Estimates in accordance with Section 4 subsection 2 showed that a fairly wide band of the spectrum of the longitudinal oscillations wsll from n L 10 to 20) and a narrower band of radial oscillations (w,, from n 5 4 to 8) were resolved without substantial errors in the calculations carried out. 2. In the first series of calculations (Re, = 102, 103, 104) in the statement of the boundary conditions on rl, r2 the model {m, p} = {,%, (0)) S,(L) } was used: pu, (0, r, t) =0.2 (l-r’) T(0,
r, t) =I,
[ l-exp
(-0.2t)]
20,
~~(0, r, t) =O;
(5.1) p
CL, r, t>=p,=lly;
T(L,
r, t)=1,
uz(L,
r, t) =0
if u,CO on r2,
rE [0, I].
155
Unsteady subsonic flows of viscousgas
As a result of the specification on I’1 of a monotonic growth of the flow rate the general development of the flow in the duct was accompanied by the generation of intensive longitudinal anharmonic oscillations with a clearly defined dominant frequency ox. 1 = 0.186, close to the fun&mental characteristic frequency of the one-dimensional quarter-wave resonator w,lo = n/2L = 0.196. Also resolved in the calculations were the radial oscillations (that is, normal waves of higher for orders [ 1,6]) of considerably less intensity ( max 1p (x, 1, t) -p (5, 0, t) 1/p,G10s3 r Z IO) with fundamental frequency w,.r = 3.6 % wxl. We note that the generation of radial oscillations arises not only from inhomogeneous conditions on rl (if dpu, (0, r, t) /drSO, then dp (0, r, t) /6’r+O), but also due to refraction effects in the vortex flow [ 11. For Re, = 102 the attenuation of the oscillations and the establishment of a steady state of the flow with p (0, r, r)/p, = 1.09 was observed for r 2 100. In all the cross sections the steady state profiles of ~1, pul were close to parabolas.
FIG. 1 ,
FIG. 2
Figure 1 shows graphs of m, (t ) /mar m2(t)lmo (m,=m,(~)=m~(m) =x./10). igs l-3 respectively) for Re, = 103, t < 60. We note that the curves of ~(0, 0, G/pa P p (0,0, t)/pa in Figs. 1-3 are averaged over small radial oscillations. For Re, = 102 in the given problem ox,,= (zlL.)=,,, (a,) X,,B (a,) Iti, at least for wxn with n 6 10. As Re, -+ 00(II + 0) for all frequencies wxn we must have low + oVI + 0, Q + or (see (4.6) (4.7)). Considering in G 2 homogeneous flow of an ideal gas with mean velocity (U I )a = MU z 0.1 and putting Qr = --au, Q2 = =(see (3.1) (3.5)) it is easy to show that or = 0 (@,/L) = 0 (10-z). Here the value of or, to a first approximation independent of frequency, is determined mainly by the convective transfer of acoustic energy through the cross section r2, From the numerical solution of the problem for the fundamental tone wX1 it is possible to determine the effective attenuation factors (YZ 0.6, a = 0.03, (Ya 0.02 (for Re, = 102, 103, 104 respectively), which agree with the estimates made. For Re, = 104 (Re, = RecMmax = 2000) an unsteady separation of the flow from the side wall, the development of secondary processes on instability in the viscous shear layer being torn off etc. were observed. These precritical phenomena are discussed individually in [2,3].
156
A. T. Fedorchenko
3. Also solved was a version of the problem (for Re, = 103) in which only one change in (5.1) was made: instead of specifying on I’1 the temperature distribution T(0, r, t) = 1 a condition of the type (2.6) is introduced: T(0,
whereg(r)=
r, t)=l+g(r)
1 forrE
[ (~(0,
[0, l-S],g(r)=
r. t)lp,)(i-i)‘Y-l],
1 - [(r-
rE[O,
1 t&)/S]2 forrE
11,
t4,
[l -6, 11.
This condition is equivalent to specifying a quasi-stationary entropy profile on the segment rl: (x = 0,r E [0,1 -61). In the calculations 6 = 0.2 was used (that is, 6=0 (6,“‘) =0 (a,‘“‘), see (4.7)). As is to be expected the amplitudes of p (0, 0, t)/~~ exceeded only slightly (by not more than 10%) the corresponding values of the version (5.1), the general picture of the course of the process being practically similar. 4. In another version, Re, = 103, in relation to (5.1) only the condition on r2 was changed in accordance with the model (3.4): P c-k I’,
where@(t)=OfortE
t) =Pl+p (t) [m(t) --m* (t) 1,
rE[O,
I],
t4,
[0,tb],P(t)~00.28fort>tb.
The instant ~b = 11.7 corresponds to the first equalization of the outflows ml and m2 in Fig. 1. Therefore, by (3.9, (3.6), on r2 in the linear approximation we have s22 = 00(92 = -1) for t < &, Q2 = 1 (172 14 1) for t > tb, that is, for c > cb there is practically no reflection of the quasiplane wave from the cross section r2, which must lead to the rapid attenuation of the longitudinal oscillations in G. Figure 2 for t Q 60 shows graphs of m 1 (t)/ma, m2 (t)/ma, p (0, 0, t)/pa, p (15, 0, t)/po (curves l-4 respectively). It is seen that for the solution of the non-linear problem the results of the linear estimates for Ql, Q2 are confirmed. The attenuation of the oscillations and the establishment of the steady state of the flow is observed in practice for f 2 50. It is obvious that such an effective attenuation of the longitudkl oscillations (in the entire frequency spectrum computed) could not have been accomplished by introducing some model of volume absorption. 5. A similar attenuating effect was obtained on solving the problem with Re, = 103 also, where unlike (5.1), on rl a model S,, of the form (2.7), (2.8) was used: Pf
(0, r, t>=cp(r) be”(t) --l-f CL)1 t
ph(O,
r, t) =O,
T(0,
r, t)=l
if
u1 (0, r, t) 20,
(5.2) peo(t)=pa{l+0.2[1-exp
(-0.2t)]},
cp(r)=0.2(1-rZ)/0.191p,~=y(1-r’),
p”(t)=p(O,
0, t),
rE[O,
t=Q.
I],
From (5.2). (3.3) we can estimate Qr = 0 (-I), In1 I < 1.~2 * - 1. that is. in this case also rapid attenuation of all the frequencies w,, can be expected. Figure 3 shows curves 1-4 corresponding to m, (t) lm,, m2 (t) lm,, p (0, 0, t) lp,,, Establishment of the steady state (with the outflow m, = ml = m2 differing p,“(t) /pa. negligibly for the parameters of (5.2) from m, in problem (5.1)) is observed for t>,50.
1.57
Unsteady subsonic flows of viscous gas
'a
20
4D
t
&!7
FIG. 3
6. Using the model {p, I)} the problem of the acceleration of a subsonic flow of an originally quiescent gas when there is a monotonic increase of the pressure on rl was solved numerically: ~(0, r. t)=p,{l+d[l-esp
(--h-t)]}.
T(0,
r, t)=l,
~~(0. r, t) =O, p(L,
ford>O,k>O;ul
r, t) =pa=const
>,Oonrl,r2.
We recall that the solution of the linear problem of the acceleration of the flow of a viscous incompressible fluid in an infinitely long circular tube (on the instantaneous application of a specified longitudinal pressure gradient to the quiescent fluid) was first given in [ 131, and later (and apparently independently) in [ 141. Under real conditions the accelerated flow of a compressible gas in a duct of finite length is usually realized by varying the pressure at one end of the duct, which may cause the development of longitudinal acoustic oscillations of considerable amplitude. Here the solution of the problem for Re, = 100, d = 0.1, k = 0.2 is shown. The choice of such a low Re, number was determined by the fact that for these values of Re,, L the effects of viscosity show themselves clearly even in the course of a comparatively short time interval. The calculation of this version was continued up to ff = Re, = 100, when in practice the attenuation of the oscillations and the establishment of the steady state of the flow was observed. If we introduce the relative time t* = t/Re,, independent of the choice of the characteristic velocity co, then the total “acceleration time” tf* x 1 agrees with the solution of [ 141.
Figure 4 shows graphs of w,(t), u:,(i) (w,(t) =mj(t)/n), ~(0, T’, t)/pa, p(L/& 0, t)/p, (curves l-4 respectively). For comparison curve 5 gives the graph of wf(r*)wl (lOO)/wf (I), where wf (t*) is the outflow from the solution of [ 141. If we estimate the “switching in” time of the constant longitudinal pressure gradient of value tR * I/k * 5. then by a shift of curve 5 to the right by At * 5 it will be close to the variation of the “average” outflow (WI + ~2 )/I!.
From the curves of Fig. 4 it is possible to determine approximately for the fundamental frequency ox 1 * n/L the characteristic attenuation time tl =Z15 to 20, and also the decrease of the phase velocity (cl = 0.9), which agrees with the estimates by formulas (4.6) (4.7) (cL+O, (a,) ix l/t,=1/16> (a,) ,).
158
A. T. Fedorchenko
1
P/Pa
FIG. 4
FIG. 5 Figure 5 shows the profiles of the longitudinal velocity ~1 (0, r, t), u1 (L, r, t) at various instants of time. It is obvious that the effect of friction at the wall, as in the solutions of [ 141, is gradually transmitted along the radius to the axis. For t = 100 the profiles of ~1, pul in all cross sections are almost parabolic in shape, WI = w2 = 0.112, I (WI - wz)/wl I S 0.3%. We note that in this problem in conditions of quasi-parallel flow ( I~U~lic~]]~Zllo, (IdplLk& >~jQddyIIG) the formation of the profiles UI (r) (including that on I’l) is primarily determined by the terms apldx, (r Re,) -9 ( prau,ldr) ldr in the equation of motion along the Ox axis. It is characteristic that the asymptotic solutions (as t + -) of all the problems of subsections 2-6 define on a finite segment of the cylindrical duct some analogs of the Poiseuille flow for the case of a compressible medium. Therefore, using the models presented we can investigate the nonlinear effects of the action of acoustic oscillations (for example, in the spectrum of the characteristic resonance frequencies w xn, car,, of the domain G) on this important type of flow. The author thanks A. A. Samarskii, K. I. Artamonov and 0. S. Ryzhov for discussing the results and for valuable comments TramWed
by
J. Berry
159
Unsteady subsonic flows of viscousgas REFERENCES 1. NAYFEH, A. H., KAISER, J. E. and TELIONIS, D. P. Acoustics of aircraft engine-duct systems AIAA J., 13,2, 130-153,1975. 2. FEDORCHENKO, A. T. On non-linear acoustic effects in the flow of a viscous gas in a cylindrical duct of finite length. Do&LA&ad Nauk SSSR, 237,6, 1315-1318, 1977. 3. FEDORCHENKO, A. T. On calculations of two-dimensional unsteady flows of a viscous gas in a short cylindrical dust with an end blast. Izrl. Akad. Nauk SSSR. Mekhan. zhtdkosti igaza, No. 1,9-17,
1919.
4. THOMAS, P. D. Boundary conditions for implicit solutions to the compressible Navier-Stokes equations in finite computational domains. AM4 Comput. Fluid Dynamics Conj: Williamsburg Vu, 1979, pp. 14-26. 5. ROZHDENSTVENSKII, B. L. and YANENKO, N. N. Systems of quasilinearequations and their applications to gas dynamics (Sistemy kvazilineinykh uravnenli i ikh prilozheniya k gazovoi dinamike), Nauka, Moscow, 1978. 6. ISAKOVICH, M. A. Generalacoustics (Obshchaya akustika), Nauka, Moscow, 1973. 7. LANDAU, L. D. and LIFSHITS, E. M. Mechanics of continuous media (Mekhanika sploshnykh sed), Gostekhizdat, Moscow, 1953. 8. FEDORCHENKO, A. T. Numerical investigation of some MHD flows of a viscous compressible gas in a rectangular cavity. Izv. A&ad. Nauk SSSR. Mekhan. zhidkosti igaza, No. 5,27-33, 1975. 9. FEDORCHENKO, A. T. Numerical investigation of an acoustic resonance in a rectangular cavity in a viscous gas flow. Akustich zh., 24,5, 75 l-759, 1978. 10. SAMARSKII, A. A. and POPOV, Yu P. The difference schemes of gas dynamics (Raznostnye dinamikil. Nauka, Moscow, 1975. 11. ROACHE, P. J. Computationalfluid dynamics Hermosa PubL, Albuquerque,
skhemy gazovoi
1976.
12. RAYLEIGH. Theory of sound (Teoriya zvuka), VoL 2, Izd-vo in lit, Moscow, 1955. 13. GROMEKA, I. S. On the theory of motion of a f7uid in narrow cylindrical tubes (K teorii dvizheniya zhidkosti v uzkikh tsilindricheskikh trubkakh). Kazan, 1882. 14. SZYMANSKY. P. Quelques solutions exactes des equations de I’hydrodynamique de fluide visqueux dans le cas d’un tube cylindrique. J. Math. pures et appl, ser. 9, 11, 1, 67-107, 1932.