- Email: [email protected]
TVD schemes for the calculation of flow in pipes of variable cross-section
Mathl. Comput. Modelling Vol. 21, No. 3, pp. 85-92, 1995
Pergamon
08957177(94)00216-9
Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
TVD Schemes for the Calculation of Flow in Pipes of Variable Cross-Section Department
Department
J. M. CORBERAN of Applied Thermodynamics, Universidad 46071 Valencia, Spain M. LL. GAS&N of Applied Mathematics, Universidad 46071 Valencia, Spain
Politkcnica
Politecnica
(Received June 1994; accepted September
de Valencia
de Valencia
1994)
Abstract-In
this paper, first and second order explicit TVD finite difference schemes have been adapted to the calculation of the unsteady one-dimensional flow in ducts of varying cross-sectional area. The basic idea to bear in mind has been the formulation of the schemes in function of the flux instead of the conserved variables, in order to include the term concerning the cross-section variation correctly. The calculation of eigenvalues of Jacobian matrix has been baaed on the Roe’s linearization technique, and a modification of the first and third eigenvalues has been done in order to force the satisfaction of the entropy condition. Numerical results are presented to demonstrate the performance of the proposed scheme. Keywords-Conservation
laws, TVD scheme, Entropy condition,
Riemann problem
1. INTRODUCTION To find good approximations of the discontinuous solutions been a constant objective in modern numerical mathematics.
of systems of conservation laws has On the one hand, upwind schemes
have been used over the years because of their success solving this problem [l], but they have a strong numerical dissipation which spreads discontinuities over many grid points, and it produces a low accuracy in the smooth region of the solution. On the other hand, it is well known that the classical second-order methods, such as Lax-Wendroff, give good resolution in smooth regions but develop an excessive oscillation near shock waves and steep gradients [2]. To avoid the above difficulty, some new shock capturing finite difference schemes have been constructed, these methods are of second-order accuracy in the smooth region of the solution without presenting the spurious oscillations associated with the conventional second-order schemes in the presence of discontinuities. An important kind of these schemes was introduced by Harten [3]; they are the TVD (total variation diminishing) schemes, which guarantee that monotone profiles remain monotone. He derived conditions on the coefficients for explicit and implicit fully disthe TVD property and second-order accuracy for a scalar crete schemes [3-51, which guarantee approximation. After having designed a TVD scheme in the scalar case, he extended it formally to systems, using a field-by-field limiter and Roe’s decomposition. Since then many other second-order accuracy TVD schemes have been proposed [6,7]. The authors would like to thank the DGICYT support to the realization of this work.
(Direccibn
General de Investigacibn
CientSca
y Tecnica)
‘b-et 85
for their
by AH-W
J. M. CORBER~
86
AND M. LL. GASC~N
Special attention should be made to the proper formulation of the above mentioned schemes when they are applied to particular cases. In this paper, we describe the development of TVD schemes to calculate one-dimensional flow in pipes of variable cross-section. In this case, not only it is sufficient to add an integration of the term corresponding to the cross-section variation, but it is necessary to do an especial adaptation of the TVD technique to this concrete problem.
The basic idea to bear in mind is the formulation of the schemes with a definition of
fluxes including the term concerning the cross-section variation, in order to secure conservation. The organization of this paper is the following: In Section 2, we describe the system of conservation laws that governs the one-dimensional unsteady flow in pipes of variable cross-section, and a new integration for the term concerning the cross-section is presented.
In Section 3, an
adequate formulation of a first order TVD scheme is given. In Section 4, the development of a second-order
scheme, for regions where the solution is smooth, is presented.
This scheme is a
good compromise to solve this type of flow because it leads to good results and allows to take into account the effect of the cross-section variation without the use of 2-D or 3-D schemes, which should be of course much more time consuming. It is consistent with the Roe’s scheme in the constant diameter case, and therefore, it may admit stationary nonphysical discontinuities, so a modification of this scheme, based on the choice of the first and third eigenvalues of Jacobian matrix, has been introduced in the new scheme without disturbing the conservation; it allows the calculation of the correct physical solutions, i.e., solutions satisfying the entropy condition [8]. Finally, numerical experiments using the previous schemes are presented in the last section.
2. EQUATIONS The partial differential equations governing one dimensional unsteady nonhomentropic pipes of variable cross-section are written as: Wt + F, + B’ = C,
flow in
(1)
where
(2) and
(3) with u = velocity, p = pressure, p = density, A = cross-section, y = ratio of specific heats, q = term of heat transferred to the walls, g = friction term, t = time, 2 = distance, Al corresponds to the wall surface per unit of length, so that corresponds to the calculated one by the hydraulic diameter. The Jacobian matrix J(W) of the above system is independent of the cross-section and is given by 0 J(W)
=
1
$7.9 i +(y - l)G
-47 -uH
H-
0 - 3)
(y-l)U2
Y-l yu
7
(4
1
where H is the enthalpy, with pH = E Sp and E = pu2/2 +p/(y - 1). The flux terms are the natural fluxes that correspond to mass and energy in 1-D flow, whereas the momentum equation establishes the equilibrium among the different forces that act on the control volume.
TVD Schemes
87
The pdA term represents the reaction forces from the wall to the flow in pipes where the section of the product pdA over the control volume:
varies and its evaluation requires the integration Bc2! = _ 293
J
Note that this integration
i
’ pdA = (b’pi + (1 - 0)~~) (Ai - Aj) ,
8 E [O,11.
consists of a linear function of pi and pj with constant
(5) weighting
factors; several definitions for the parameter 0 have been studied, although an average of pressure, i.e., 6’ = l/2 (trapezoidal
rule integration),
produces good results.
Now, Bi,j denotes the following vector Bij=(O
Bc2! O)T %>3
3. THE TVD
(6)
SCHEME
In this section, we present a TVD scheme for the system of conservation laws (l), including the integration of the term concerning the cross-section variation. Firstly, the TVD scheme is defined for the scalar equation and then we extend the scheme to the whole system of conservation laws. For the scalar conservation law: Wt i f(w),
= 0
(7)
a 3-point finite difference scheme in conservative form becomes given by w;+l
= wa - x fin+l,2 [
fin-l/2I I
(8)
where the flux in the middle could be defined as
fin+l/Z=
’ [f:+fir;1 -
2
with
x 4+1/2
=
(G+I,~) (G+I - w:)]
;Q
f%il-f% w,+1-w,
ai,
’
wz+1 -
wi # 0,
Wifl
W, =
-
0.
(10)
Harten showed in [3] that the above scheme is TVD if Qsatisfies the inequalities In] < Q(Q) < 1, under the CFL restriction /CX] 5 1. In the variable cross-section case, it will be necessary to avoid the differences on w to introduce the term concerning the cross-section variation in the middle flux, therefore we have modified the half flux in the TVD scheme for the scalar equation as:
(11) and it is possible to show that the scheme is TVD under the CFL-restriction the following inequalities: 1 < h(a) < 1 ’ 0 1 (1.-< h(a) 5 -1,
O
Q < 0.
The choice of the function h for our upstream scheme has been
which assures the monotonicity
property.
(~1 5 1, if h satisfies
(12)
J. M. CORBEF&N AND M. LL. GA&N
88
Let us start with the extension of this TVD scalar scheme to the system of conservation laws (1). The Jacobian matrix J(W) is diagonalizable, the diagonal matrix A is the eigenvalues matrix, whose main diagonal is defined by u - a, u and u + a, respectively; the associated right eigenvectors matrix (R) is 1 R=
1
u-a ( H - ua
1
u u+a , $2~~ H +ua 1
(14)
where a is the sound speed, and the left eigenvectors matrix (L) is
(
&++a)
L=
5
)
l(f$--u2)
( If J(W)
-(s
ua
4
(
2uatua
a -2a
2 -c - ua 7-i > (
iu - a 7-l )
.
(15)
a 1
E J, i.e., independent of z and t, equation (1) is the same as LWt + ALW, + LB’ = LC
(1’3)
and if we change the variable W by LW then it is possible to analyse every equation separately. This method is W?+r=W~?-X F.+ -F: a-1/2 + (At)C$ (17) Z [ 2+1/2 I where the vector flux is given by: FiL,2
= f
[cFF
+ Fin;11 *
(Bti+l/2
+ B,n+l/2,i+l) (18)
-
(Rh(h)L)Y+l/2
(
F,“+1 - Fz? + BTi+l/2
+ %+1/2,i+l
)
+ At(RAL)r+l/2CF+1/2
1
with h(A)i+1/2, a diagonal matrix of which the main diagonal is defined by sign(u - u)~+~,~, sign(u)i+l,z and sign(u+o)i+r/2, respectively. Note that the evaluation of the Bi,j terms requires a double estimation of the pressure at the grid points and in the middle of the grid points, therefore the pressure field must have a double definition. The calculation of eigenvalues in the middle is based on the existence of a mean value Jacobian Ji+1/2 z J(Wi Wi+l), based on the Roe’s linearization [9], which satisfies 1
Fi+i - Fi = J (Wi, Wi+l) (Wi+l J (Wi, Wi) = Ji,
J (Wi+l, Wi+l)
- Wi) (19)
= Ji+l*
For our system of equations, this averaging takes the following form:
x+1/2
%+1/2 =
Xi+l/ZUi+l
+ Ui
Xi+1/2 + l
%+I/2
=
>
f&+1/2
=
Xi+l/2Hi+l+ Xi+1/2 + l
Hi
(20)
=
The previous scheme is a total variation diminishing (TVD) scheme under the CFL restriction,
X max { J%+1/2\ + i
%+1/2}
5 1,
(21)
which guarantees that the scheme does not generate spurious oscillations. It is easy to see that the flux defined by (18) is similar to the Roe-flux in the constant section case, and this is a good choice for the flow calculation.
TVD Schemes
4. EXTENSION
TO SECOND
89
ORDER
ACCURACY
In this section, we deal with the extension to second order of the scheme described previously. Our extension technique is based on the procedure suggested by Harten 131, but for our case, we express every difference on w as flux differences. The modified equation of the scalar scheme is w + fz = & (la] - a2) (Az>w,,
+ g
(31~4 - 1 - 2cu2) (AE)~w~~~ + 0 (Az3)
(22)
and, for the linear case, can be transformed on wt + f3c = ;(h(u)
- a)(Az)fzz
+ $ (3101 - 1 - 2cu2) (A~)~fi,,
where
I,
h(cr) = Therefore,
+ 0 (AZ”) ,
Q L 0,
-1,
(23)
(24
cr
the modified equation can be rewritten as w +
f - +(4
= ; (31~4 - 1 - 2a2) (Ax:)~~~~, + 0 (Ax’).
- WWfz]
z The argument to transform the previous scheme in a second-order approximation of the above first order scheme to wt + [f + Then, to approximate
(25)
scheme is based on the
(26)
41% = 0.
the flux c$, we define 6+1/2
= ;
(h (4+1,2)
-
(4+1/z))
(fi71+1
-
(27)
fi”)
and
with sY+~,~ = sign[$y++,,2]. Now, if 4% l-4.
pz”+l,2
=
{
o”.:l-w.’ wi+l Wi+l
7
-
wi+O7 Wi
(29)
O
=
and (30) Then, the new second order flux for the scalar case is defined by fiTbl/2 = 5’ [f,” + fi’“+l +
4; + 4:+“,, - h (a:+,,,) (.fi7;1 - fi^>- h” @+1,2)
(4:+;1- 431 ) (31)
which can be reduced to
and finally, the flux is given by fGl12
= ;
[f;
+ f;+l
-
h @+1,2)
(&
-
f:)]
+ (Py+1/23
(33)
J. M. CORBEF&N AND M. LL. GAS&N
90
where (34)
y$“+rlz = ~~+i,~ max { Therefore,
the formulation of the new TVD flux of second order for the calculation can be written as
of flows in
pipes of variable cross-section
F&l,2
=
(Fin + F$)
f
f
( q+1/2
-
+ Bzn+1/2,i+1
ww-x+1/2
+
(
where
4?+1/2
=
;
(h
(K+1,2)
-t %/2,i+1
1
(35)
)
+ %‘+1/24’+1/27
and A are matrices described
R, L,h(A)
in
of the vector
( w:+1,2))
-
Fl? + qi+1/2
of the second-order flux vector is
denotes the j-component
($i+l/2)j
JTi-1 -
AWW++,,2CF+l/2
where F, B and C are vectors defined previously; Section 3 and the j-component
)
LY+1/2
(FiFl
-
Ft? + qi+1/2
+ B,;1,2,i+1)
*
(37)
This scheme allows the solution of shocks with high resolution, but it also may admit entropy violating discontinuities as solution. In the following, a modification of this scheme is given to force the satisfaction
of the entropy condition.
The entropy inequalities
introduced by Lax [8] are Ui+l
<
V <
Ui+i + aifi
<
V < Ui + Ui,
Ui+l
-
Ui
-
Ui,
(38)
respectively, for the first and third eigenvalues of the Jacobian matrix, where v is the shock speed. In order to secure the previous condition, we may modify the numerical signal speeds in the average of the Jacobian matrix by V!%+1/2 =
min{%+l/2
v&1,2
=
U&l/2
V&l,2
=
max
-
%+1/27
ui -
%) (39)
{%+1/2
+ ai+l/2,
Ui+l
.
+ %+l}
And the matrix h(A) is a diagonal matrix of which the diagonal’s elements are
h
(
&l/2
h (Vf+l,2)
)I =
=
sign(~uzf+~,~ >, if vZ~+,,~=
sign
‘Ux-az
~.+1/2-%+1/2 (VZ:l,2)
T
if Vuz!+1,2
Ui+1/2
-
%+1/2,
=ui-ai,
(40)
Therefore, the numerical flux can be calculated by the equation (35) where the h(A) matrix has been replaced by the above expression, resulting an entropy-satisfying scheme.
TVD Schemes
91
Figure 1. Pipe with 1 m length, 0.05 m diameters at both sides and 0.038 m diameter at the throat.
0.5
-
I
0 0
0.2
0.4
0.6
Figure 2. Pressure solution in the convergent-divergent
5. NUMERICAL In this section,
1
08
nozzle.
RESULTS
we show first the numerical results of pressure obtained
with the proposed
second-order TVD scheme, for the calculation of the flow in a convergent-divergent ure 1 shows the geometrical dimensions of the nozzle.
nozzle. Fig-
The calculations have been performed with 2 bar of pressure at the left side and 1.5 bar at the right side and 300 K of temperature at both sides. The flow at the entry is supposed homentropic whereas at the exit the pressure is that of the atmosphere. Figure 2 shows the numerical results of pressure for the convergent-divergent nozzle when the steady solution is carried out, obtained with a 42 uniform grid and under the CFL restriction of 0.95. The exact solution in Figure 2 is shown by the solid line and the values of the numerical solution are indicated by symbols. The numerical solution shows very good conservation, almost negligible diffusion and gives a reasonable resolution of the shock even for a lower number of nodes. Figure 3 shows the numerical solution of the pressure after 3.5 seconds for the Riemann problem with a nonphysical discontinuity. This problem, described in 131, is given by the following initial conditions (41)
HCM 21:3-G
J. M. COFI.EIEFI,.&N AND M. LL. GAS&N
92 20
I
n n n n n n
m I
Figure 3. Pressure solution obtained
for the Riemann problem (41).
The calculations in Figure 3 were performed with 140 cells and the same CFL restriction as the above problem. As can be observed, the scheme with the modification (40) gives an entropysatisfying solution. The authors have tried to adapt the FCT techniques (Flux-Corrected Transport) [lo] and the original TVD scheme proposed by Harten [3] to the calculation of cases with variable crosssection. The most serious difficulty has been found in the treatment of the differences expressed on the variation of the conserved variables W. The terms concerning those W differences led always to excessive diffusion in the points in which the rate of cross-section variation is high, or even to problems in the global conservation along the pipe. This difficulty disappears when all the differences in the scheme can be expressed as flux differences. The proposed scheme has been derived following this idea, but this has been only possible for a TVD scheme and not for the FCT techniques. The second-order TVD finite difference scheme explained in this paper has been used to calculate the flow in pipes of Internal Combustion Engines, having shown always a very good behaviour in the presence of the thermal and pressure discontinuities characteristic of this kind of engine.
REFERENCES 1. B.V. Leer, On the relation between the upwind-differencing schemes of Godunov, Engquist Osher and Roe, Siam J. Sci. Stat. Comput. 5 (1) (1984). 2. A. Sod, Numerical Methods in Fluid Dynamics, Cambridge University Press, (1985). 3. A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49, 357-393 (1983). finite-difference schemes, Siam J. Numer. 4. A. Harten, On a class of high resolution total-variation-stable Anal. 21, l-23 (1984). 5. H.C. Yee, R.F. Warming and A. Harten, Implicit total variation (TVD) schemes for steady-state calculations, J. Comput. Phys. 57, 327-360 (1985). 6. P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, Siam J. Numer. Anal. 21, 995-1011 (1984). 7. P.K. Sweby, High resolution schemes using flux limiters, Lectures in Applied Mathematics, Volume 22, pp. 289-309, (1985). a. P.D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, SIAM, Philadelphia, (1972). 9. P.L. Hoe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43, 357-372 (1981). 10. D.L. Book, J.P. Boris and K. Hain, Flux-corrected transport II: Generalizations of the method, J. Comp. Phys. 18, 248-283 (1975).
Pergamon
08957177(94)00216-9
Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
TVD Schemes for the Calculation of Flow in Pipes of Variable Cross-Section Department
Department
J. M. CORBERAN of Applied Thermodynamics, Universidad 46071 Valencia, Spain M. LL. GAS&N of Applied Mathematics, Universidad 46071 Valencia, Spain
Politkcnica
Politecnica
(Received June 1994; accepted September
de Valencia
de Valencia
1994)
Abstract-In
this paper, first and second order explicit TVD finite difference schemes have been adapted to the calculation of the unsteady one-dimensional flow in ducts of varying cross-sectional area. The basic idea to bear in mind has been the formulation of the schemes in function of the flux instead of the conserved variables, in order to include the term concerning the cross-section variation correctly. The calculation of eigenvalues of Jacobian matrix has been baaed on the Roe’s linearization technique, and a modification of the first and third eigenvalues has been done in order to force the satisfaction of the entropy condition. Numerical results are presented to demonstrate the performance of the proposed scheme. Keywords-Conservation
laws, TVD scheme, Entropy condition,
Riemann problem
1. INTRODUCTION To find good approximations of the discontinuous solutions been a constant objective in modern numerical mathematics.
of systems of conservation laws has On the one hand, upwind schemes
have been used over the years because of their success solving this problem [l], but they have a strong numerical dissipation which spreads discontinuities over many grid points, and it produces a low accuracy in the smooth region of the solution. On the other hand, it is well known that the classical second-order methods, such as Lax-Wendroff, give good resolution in smooth regions but develop an excessive oscillation near shock waves and steep gradients [2]. To avoid the above difficulty, some new shock capturing finite difference schemes have been constructed, these methods are of second-order accuracy in the smooth region of the solution without presenting the spurious oscillations associated with the conventional second-order schemes in the presence of discontinuities. An important kind of these schemes was introduced by Harten [3]; they are the TVD (total variation diminishing) schemes, which guarantee that monotone profiles remain monotone. He derived conditions on the coefficients for explicit and implicit fully disthe TVD property and second-order accuracy for a scalar crete schemes [3-51, which guarantee approximation. After having designed a TVD scheme in the scalar case, he extended it formally to systems, using a field-by-field limiter and Roe’s decomposition. Since then many other second-order accuracy TVD schemes have been proposed [6,7]. The authors would like to thank the DGICYT support to the realization of this work.
(Direccibn
General de Investigacibn
CientSca
y Tecnica)
‘b-et 85
for their
by AH-W
J. M. CORBER~
86
AND M. LL. GASC~N
Special attention should be made to the proper formulation of the above mentioned schemes when they are applied to particular cases. In this paper, we describe the development of TVD schemes to calculate one-dimensional flow in pipes of variable cross-section. In this case, not only it is sufficient to add an integration of the term corresponding to the cross-section variation, but it is necessary to do an especial adaptation of the TVD technique to this concrete problem.
The basic idea to bear in mind is the formulation of the schemes with a definition of
fluxes including the term concerning the cross-section variation, in order to secure conservation. The organization of this paper is the following: In Section 2, we describe the system of conservation laws that governs the one-dimensional unsteady flow in pipes of variable cross-section, and a new integration for the term concerning the cross-section is presented.
In Section 3, an
adequate formulation of a first order TVD scheme is given. In Section 4, the development of a second-order
scheme, for regions where the solution is smooth, is presented.
This scheme is a
good compromise to solve this type of flow because it leads to good results and allows to take into account the effect of the cross-section variation without the use of 2-D or 3-D schemes, which should be of course much more time consuming. It is consistent with the Roe’s scheme in the constant diameter case, and therefore, it may admit stationary nonphysical discontinuities, so a modification of this scheme, based on the choice of the first and third eigenvalues of Jacobian matrix, has been introduced in the new scheme without disturbing the conservation; it allows the calculation of the correct physical solutions, i.e., solutions satisfying the entropy condition [8]. Finally, numerical experiments using the previous schemes are presented in the last section.
2. EQUATIONS The partial differential equations governing one dimensional unsteady nonhomentropic pipes of variable cross-section are written as: Wt + F, + B’ = C,
flow in
(1)
where
(2) and
(3) with u = velocity, p = pressure, p = density, A = cross-section, y = ratio of specific heats, q = term of heat transferred to the walls, g = friction term, t = time, 2 = distance, Al corresponds to the wall surface per unit of length, so that corresponds to the calculated one by the hydraulic diameter. The Jacobian matrix J(W) of the above system is independent of the cross-section and is given by 0 J(W)
=
1
$7.9 i +(y - l)G
-47 -uH
H-
0 - 3)
(y-l)U2
Y-l yu
7
(4
1
where H is the enthalpy, with pH = E Sp and E = pu2/2 +p/(y - 1). The flux terms are the natural fluxes that correspond to mass and energy in 1-D flow, whereas the momentum equation establishes the equilibrium among the different forces that act on the control volume.
TVD Schemes
87
The pdA term represents the reaction forces from the wall to the flow in pipes where the section of the product pdA over the control volume:
varies and its evaluation requires the integration Bc2! = _ 293
J
Note that this integration
i
’ pdA = (b’pi + (1 - 0)~~) (Ai - Aj) ,
8 E [O,11.
consists of a linear function of pi and pj with constant
(5) weighting
factors; several definitions for the parameter 0 have been studied, although an average of pressure, i.e., 6’ = l/2 (trapezoidal
rule integration),
produces good results.
Now, Bi,j denotes the following vector Bij=(O
Bc2! O)T %>3
3. THE TVD
(6)
SCHEME
In this section, we present a TVD scheme for the system of conservation laws (l), including the integration of the term concerning the cross-section variation. Firstly, the TVD scheme is defined for the scalar equation and then we extend the scheme to the whole system of conservation laws. For the scalar conservation law: Wt i f(w),
= 0
(7)
a 3-point finite difference scheme in conservative form becomes given by w;+l
= wa - x fin+l,2 [
fin-l/2I I
(8)
where the flux in the middle could be defined as
fin+l/Z=
’ [f:+fir;1 -
2
with
x 4+1/2
=
(G+I,~) (G+I - w:)]
;Q
f%il-f% w,+1-w,
ai,
’
wz+1 -
wi # 0,
Wifl
W, =
-
0.
(10)
Harten showed in [3] that the above scheme is TVD if Qsatisfies the inequalities In] < Q(Q) < 1, under the CFL restriction /CX] 5 1. In the variable cross-section case, it will be necessary to avoid the differences on w to introduce the term concerning the cross-section variation in the middle flux, therefore we have modified the half flux in the TVD scheme for the scalar equation as:
(11) and it is possible to show that the scheme is TVD under the CFL-restriction the following inequalities: 1 < h(a) < 1 ’ 0 1 (1.-< h(a) 5 -1,
O
Q < 0.
The choice of the function h for our upstream scheme has been
which assures the monotonicity
property.
(~1 5 1, if h satisfies
(12)
J. M. CORBEF&N AND M. LL. GA&N
88
Let us start with the extension of this TVD scalar scheme to the system of conservation laws (1). The Jacobian matrix J(W) is diagonalizable, the diagonal matrix A is the eigenvalues matrix, whose main diagonal is defined by u - a, u and u + a, respectively; the associated right eigenvectors matrix (R) is 1 R=
1
u-a ( H - ua
1
u u+a , $2~~ H +ua 1
(14)
where a is the sound speed, and the left eigenvectors matrix (L) is
(
&++a)
L=
5
)
l(f$--u2)
( If J(W)
-(s
ua
4
(
2uatua
a -2a
2 -c - ua 7-i > (
iu - a 7-l )
.
(15)
a 1
E J, i.e., independent of z and t, equation (1) is the same as LWt + ALW, + LB’ = LC
(1’3)
and if we change the variable W by LW then it is possible to analyse every equation separately. This method is W?+r=W~?-X F.+ -F: a-1/2 + (At)C$ (17) Z [ 2+1/2 I where the vector flux is given by: FiL,2
= f
[cFF
+ Fin;11 *
(Bti+l/2
+ B,n+l/2,i+l) (18)
-
(Rh(h)L)Y+l/2
(
F,“+1 - Fz? + BTi+l/2
+ %+1/2,i+l
)
+ At(RAL)r+l/2CF+1/2
1
with h(A)i+1/2, a diagonal matrix of which the main diagonal is defined by sign(u - u)~+~,~, sign(u)i+l,z and sign(u+o)i+r/2, respectively. Note that the evaluation of the Bi,j terms requires a double estimation of the pressure at the grid points and in the middle of the grid points, therefore the pressure field must have a double definition. The calculation of eigenvalues in the middle is based on the existence of a mean value Jacobian Ji+1/2 z J(Wi Wi+l), based on the Roe’s linearization [9], which satisfies 1
Fi+i - Fi = J (Wi, Wi+l) (Wi+l J (Wi, Wi) = Ji,
J (Wi+l, Wi+l)
- Wi) (19)
= Ji+l*
For our system of equations, this averaging takes the following form:
x+1/2
%+1/2 =
Xi+l/ZUi+l
+ Ui
Xi+1/2 + l
%+I/2
=
>
f&+1/2
=
Xi+l/2Hi+l+ Xi+1/2 + l
Hi
(20)
=
The previous scheme is a total variation diminishing (TVD) scheme under the CFL restriction,
X max { J%+1/2\ + i
%+1/2}
5 1,
(21)
which guarantees that the scheme does not generate spurious oscillations. It is easy to see that the flux defined by (18) is similar to the Roe-flux in the constant section case, and this is a good choice for the flow calculation.
TVD Schemes
4. EXTENSION
TO SECOND
89
ORDER
ACCURACY
In this section, we deal with the extension to second order of the scheme described previously. Our extension technique is based on the procedure suggested by Harten 131, but for our case, we express every difference on w as flux differences. The modified equation of the scalar scheme is w + fz = & (la] - a2) (Az>w,,
+ g
(31~4 - 1 - 2cu2) (AE)~w~~~ + 0 (Az3)
(22)
and, for the linear case, can be transformed on wt + f3c = ;(h(u)
- a)(Az)fzz
+ $ (3101 - 1 - 2cu2) (A~)~fi,,
where
I,
h(cr) = Therefore,
+ 0 (AZ”) ,
Q L 0,
-1,
(23)
(24
cr
the modified equation can be rewritten as w +
f - +(4
= ; (31~4 - 1 - 2a2) (Ax:)~~~~, + 0 (Ax’).
- WWfz]
z The argument to transform the previous scheme in a second-order approximation of the above first order scheme to wt + [f + Then, to approximate
(25)
scheme is based on the
(26)
41% = 0.
the flux c$, we define 6+1/2
= ;
(h (4+1,2)
-
(4+1/z))
(fi71+1
-
(27)
fi”)
and
with sY+~,~ = sign[$y++,,2]. Now, if 4% l-4.
pz”+l,2
=
{
o”.:l-w.’ wi+l Wi+l
7
-
wi+O7 Wi
(29)
O
=
and (30) Then, the new second order flux for the scalar case is defined by fiTbl/2 = 5’ [f,” + fi’“+l +
4; + 4:+“,, - h (a:+,,,) (.fi7;1 - fi^>- h” @+1,2)
(4:+;1- 431 ) (31)
which can be reduced to
and finally, the flux is given by fGl12
= ;
[f;
+ f;+l
-
h @+1,2)
(&
-
f:)]
+ (Py+1/23
(33)
J. M. CORBEF&N AND M. LL. GAS&N
90
where (34)
y$“+rlz = ~~+i,~ max { Therefore,
the formulation of the new TVD flux of second order for the calculation can be written as
of flows in
pipes of variable cross-section
F&l,2
=
(Fin + F$)
f
f
( q+1/2
-
+ Bzn+1/2,i+1
ww-x+1/2
+
(
where
4?+1/2
=
;
(h
(K+1,2)
-t %/2,i+1
1
(35)
)
+ %‘+1/24’+1/27
and A are matrices described
R, L,h(A)
in
of the vector
( w:+1,2))
-
Fl? + qi+1/2
of the second-order flux vector is
denotes the j-component
($i+l/2)j
JTi-1 -
AWW++,,2CF+l/2
where F, B and C are vectors defined previously; Section 3 and the j-component
)
LY+1/2
(FiFl
-
Ft? + qi+1/2
+ B,;1,2,i+1)
*
(37)
This scheme allows the solution of shocks with high resolution, but it also may admit entropy violating discontinuities as solution. In the following, a modification of this scheme is given to force the satisfaction
of the entropy condition.
The entropy inequalities
introduced by Lax [8] are Ui+l
<
V <
Ui+i + aifi
<
V < Ui + Ui,
Ui+l
-
Ui
-
Ui,
(38)
respectively, for the first and third eigenvalues of the Jacobian matrix, where v is the shock speed. In order to secure the previous condition, we may modify the numerical signal speeds in the average of the Jacobian matrix by V!%+1/2 =
min{%+l/2
v&1,2
=
U&l/2
V&l,2
=
max
-
%+1/27
ui -
%) (39)
{%+1/2
+ ai+l/2,
Ui+l
.
+ %+l}
And the matrix h(A) is a diagonal matrix of which the diagonal’s elements are
h
(
&l/2
h (Vf+l,2)
)I =
=
sign(~uzf+~,~ >, if vZ~+,,~=
sign
‘Ux-az
~.+1/2-%+1/2 (VZ:l,2)
T
if Vuz!+1,2
Ui+1/2
-
%+1/2,
=ui-ai,
(40)
Therefore, the numerical flux can be calculated by the equation (35) where the h(A) matrix has been replaced by the above expression, resulting an entropy-satisfying scheme.
TVD Schemes
91
Figure 1. Pipe with 1 m length, 0.05 m diameters at both sides and 0.038 m diameter at the throat.
0.5
-
I
0 0
0.2
0.4
0.6
Figure 2. Pressure solution in the convergent-divergent
5. NUMERICAL In this section,
1
08
nozzle.
RESULTS
we show first the numerical results of pressure obtained
with the proposed
second-order TVD scheme, for the calculation of the flow in a convergent-divergent ure 1 shows the geometrical dimensions of the nozzle.
nozzle. Fig-
The calculations have been performed with 2 bar of pressure at the left side and 1.5 bar at the right side and 300 K of temperature at both sides. The flow at the entry is supposed homentropic whereas at the exit the pressure is that of the atmosphere. Figure 2 shows the numerical results of pressure for the convergent-divergent nozzle when the steady solution is carried out, obtained with a 42 uniform grid and under the CFL restriction of 0.95. The exact solution in Figure 2 is shown by the solid line and the values of the numerical solution are indicated by symbols. The numerical solution shows very good conservation, almost negligible diffusion and gives a reasonable resolution of the shock even for a lower number of nodes. Figure 3 shows the numerical solution of the pressure after 3.5 seconds for the Riemann problem with a nonphysical discontinuity. This problem, described in 131, is given by the following initial conditions (41)
HCM 21:3-G
J. M. COFI.EIEFI,.&N AND M. LL. GAS&N
92 20
I
n n n n n n
m I
Figure 3. Pressure solution obtained
for the Riemann problem (41).
The calculations in Figure 3 were performed with 140 cells and the same CFL restriction as the above problem. As can be observed, the scheme with the modification (40) gives an entropysatisfying solution. The authors have tried to adapt the FCT techniques (Flux-Corrected Transport) [lo] and the original TVD scheme proposed by Harten [3] to the calculation of cases with variable crosssection. The most serious difficulty has been found in the treatment of the differences expressed on the variation of the conserved variables W. The terms concerning those W differences led always to excessive diffusion in the points in which the rate of cross-section variation is high, or even to problems in the global conservation along the pipe. This difficulty disappears when all the differences in the scheme can be expressed as flux differences. The proposed scheme has been derived following this idea, but this has been only possible for a TVD scheme and not for the FCT techniques. The second-order TVD finite difference scheme explained in this paper has been used to calculate the flow in pipes of Internal Combustion Engines, having shown always a very good behaviour in the presence of the thermal and pressure discontinuities characteristic of this kind of engine.
REFERENCES 1. B.V. Leer, On the relation between the upwind-differencing schemes of Godunov, Engquist Osher and Roe, Siam J. Sci. Stat. Comput. 5 (1) (1984). 2. A. Sod, Numerical Methods in Fluid Dynamics, Cambridge University Press, (1985). 3. A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49, 357-393 (1983). finite-difference schemes, Siam J. Numer. 4. A. Harten, On a class of high resolution total-variation-stable Anal. 21, l-23 (1984). 5. H.C. Yee, R.F. Warming and A. Harten, Implicit total variation (TVD) schemes for steady-state calculations, J. Comput. Phys. 57, 327-360 (1985). 6. P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, Siam J. Numer. Anal. 21, 995-1011 (1984). 7. P.K. Sweby, High resolution schemes using flux limiters, Lectures in Applied Mathematics, Volume 22, pp. 289-309, (1985). a. P.D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, SIAM, Philadelphia, (1972). 9. P.L. Hoe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43, 357-372 (1981). 10. D.L. Book, J.P. Boris and K. Hain, Flux-corrected transport II: Generalizations of the method, J. Comp. Phys. 18, 248-283 (1975).