Iterative Solution of Transonic Flows over Airfoils and Wings, Including Flows at Mach 1

21 Pages·2012·0.2 MB·English

Iterative Solution of Transonic Flows over Airfoils and Wings,
Including Flows at Mach 1 ;
(cid:3) y
ANTONY JAMESON
1 Introduction
Transonic aerodynamics is the focus of strong interest at the present time because it is known to encompass
one of the most e(cid:14)cient regimes of (cid:13)ight. While current transport aircraft operate just below the speed of
sound,moreradicalaerodynamicdesigns,suchasR.T.Jones’proposalforayawedwingin[1], favoroperation
at low supersonic speeds in the range from Mach 1 to 1.3 (a regime in which the sonic boom should pose
no serious problem). The hodograph method of generating a (cid:13)ow together with the corresponding boundary
shape has been perfected by Bauer, Garabedian and Korn [2]. It now provides a (cid:13)exible tool for the design of
supercriticalwingsectionswhichproduceshockfree(cid:13)owataspeci(cid:12)edspeedandangleofattack. Since,however,
it yields no information about the (cid:13)ow at o(cid:11) design conditions, and is restricted to two dimensions, it needs
to be complemented by a method for calculating the (cid:13)ow over speci(cid:12)ed shapes in two and three dimensions
throughout the desired range of speed and angle of attack. An accurate and reliable method would eliminate
the need to rely on massive wind tunnel testing, and should lead to more rational designs. Following the
initialsuccessof MurmanandCole[3], substantialprogresshasrecentlybeenmadeinthe developmentof(cid:12)nite
di(cid:11)erencemethodsforthecalculationoftransonic(cid:13)ows(cf. [4],[5],[6],[7]). Thesemethodshavegenerallybeen
restricted, however, either to (cid:13)ows which are subsonic at in(cid:12)nity, or to small disturbance theories. The present
paper describes a method using a new ‘rotated’ di(cid:11)erence scheme, which is suitable for the calculation of both
two- and three-dimensional (cid:13)ows without restriction on the speed at in(cid:12)nity, and which has been successfully
applied to a variety of (cid:13)ows over airplane wings, including the case of (cid:13)ight at Mach 1.
The mathematical di(cid:14)culties of the problem are associated primarily with the mixed hyperbolic and elliptic
type of the equations and the presence of discontinuities. The computational method should be capable of
predicting the location andstrength of the shock waves,and if it is to be capable of distinguishinga goodfrom
abadaerodynamicdesign,itshouldalsoprovideanindicationoftheassociatedwavedrag. Inthree-dimensional
applications the rapid growth in the number of points in the computational lattice also excludes any method
which is not rather economical in its use of the computer. In meeting these objectives the primary choices to
be made concern (cid:12)rst the most suitable formulation of the equations, second the construction of a favorable
coordinatesystem, andthirdthe developmentof a(cid:12)nite di(cid:11)erenceschemewhichisstable, convergent,andalso
capable of accommodating the proper discontinuities.
Since we are concerned with (cid:13)ow at fairly low supersonic Mach numbers over e(cid:14)cient aerodynamic shapes,
we may reasonably suppose that the shock waves are quite weak. A strong shock wave would cause boundary
layer separation and bu(cid:11)eting, invalidating the analysis, and a design which led to the presence of a strong
shockwavewouldinanycasebe unsuitablebecauseof the accompanyingdragriseandlossoflift. The errorin
(cid:3)TheworkforthispaperwasperformedattheCourantInstituteofMathematicalScience,NewYorkUniversity,andsupportedby
NASAunderGrantNo. 33-016-167andbytheU.S.AtomicEnergyCommissionunderContractNo. AT(11-1)-3077. Reproduction
inwholeorinpartpermittedforanypurposeoftheUnitedStatesGovernment.
yReprintedfromCOMMUNICATIONSONPUREANDAPPLIEDMATHEMATICS,VOL.XXVII,283-309(1974)
1
ignoringchangesin entropythrough shockwaves,and the resultingvorticity,should thereforebe small, andwe
canexpecttoobtainasatisfactoryapproximationbyusingthetransonicpotentialequationforirrotational(cid:13)ow
instead of the full Euler equations. In the caseof three-dimensionalcalculations, in particular,the replacement
of(cid:12)vedependentvariables(threevelocitycomponents,pressure,andenergy)byasinglevelocitypotentialleads
to important savings in machine time and memory.
Whentheboundariesareprescribed,the non-existence(cf. [8])ofsmoothtransonicsolutionsof thepotential
equation requires consideration of weak solutions (cf. [9]). The appropriate solutions admit discontinuities
acrosswhich massisconservedbut momentum isnot, sothat adragforceappears. Themethod is thus ableto
predict the onset of wave drag. If, however, we do not impose a restriction on the type of jumps to be allowed
corresponding to the entropy inequality, weak solutions of the potential equation are in general nonunique.
The equation is invariant under reversal of the (cid:13)ow direction, so that a body with fore and aft symmetry, for
example, admits both a solution with a compression shock, and a corresponding reversed (cid:13)ow solution with
an expansion shock. To obtain uniqueness we must exclude discontinuous expansions. We ought, therefore, to
restore in the numerical scheme the directional property which was removed by the exclusion of entropy from
the equations. This indicates the need to ensure that the dominant terms of the truncation error represent an
arti(cid:12)cial viscosity.
It was (cid:12)rst shown by Murman and Cole for the case of the small disturbance equations that the required
arti(cid:12)cialviscositycan be introduced in an e(cid:11)ective and simple manner by usingretardeddi(cid:11)erence formulasto
represent derivatives in the streamwise direction at all points in the hyperbolic region. This leads to a simple
implicit scheme which can be solved line by line advancing with the stream. Moreover, when the dominant
truncation error is included, the resulting di(cid:11)erential equation resembles the viscous transonic equation, which
admits solutions with the approximate structure of a shock wave (cf. [10]). The di(cid:11)erence equations exhibit
similar behavior. Shock waves emerge in the course of the calculation as compression layers spread over a few
mesh widths in which the arti(cid:12)cial viscosity becomes the dominant term.
The method described in this paper is based on a similar principle, but no assumption about the direction
of (cid:13)ow is made in constructing the di(cid:11)erence scheme. Instead, the proper directional property is obtained by
‘rotating’ the di(cid:11)erence scheme to conform with the local stream direction. This avoids the need to align one
of the coordinates with the (cid:13)ow in the supersonic zone, a source of di(cid:14)culty when this zone extends towards
in(cid:12)nityatnearsonicfreestreamMachnumbers. Therotatedschemeisderivedbyrearrangingtheequationsas
if they werewritten in astreamline and normalcoordinatesystem. Inthe hyperbolicregionretardeddi(cid:11)erence
formulas are then used for all contributions to the streamwise second derivative, producing a positive arti(cid:12)cial
viscosity correctly oriented with the stream.
The method described in this paper is based on a similar principle, but no assumption about the direction
of (cid:13)ow is made in constructing the di(cid:11)erence scheme. Instead, the proper directional property is obtained by
‘rotating’ the di(cid:11)erence scheme to conform with the local stream direction. This avoids the need to align one
of the coordinates with the (cid:13)ow in the supersonic zone, a source of di(cid:14)culty when this zone extends towards
in(cid:12)nityatnearsonicfreestreamMachnumbers. Therotatedschemeisderivedbyrearrangingtheequationsas
if they werewritten in astreamline and normalcoordinatesystem. Inthe hyperbolicregionretardeddi(cid:11)erence
formulas are then used for all contributions to the streamwise second derivative, producing a positive arti(cid:12)cial
viscosity correctly oriented with the stream.
Themethodretainsthe propertyofthe Murmanschemeofautomaticallylocatingtheshocks. Thisisagreat
advantage for the (cid:13)ows to be considered, which may contain a complex pattern of shocks above and below the
wing, together with a detached bow shock if the (cid:13)ow is supersonic at in(cid:12)nity.
2
2 Formulation of the Equations
In the isentropic (cid:13)ow of a perfect gas the equation of state
p
(2.1) =constant
(cid:26)(cid:13)
is assumed to hold, where p is the pressure, (cid:26) the density, and (cid:13) the ratio of speci(cid:12)c heats. If the density is
normalized to the value unity at in(cid:12)nity, and M. is the Mach number at in(cid:12)nity, the speed of sound a is given
by
dp (cid:26)(cid:13)(cid:0)1
(2.2) a2 = =
d(cid:26) M12
The energy equation can then be expressed as
(cid:13) 1
(2.3) a2+ (cid:0) q2 =constant
2
q being the velocity vector. Also conservation of mass requires
(2.4) :((cid:26)q)=0
r
Since the (cid:13)ow is isentropic it is also irrotational. Introducing a velocity potential (cid:30), the energy and mass
equations can then be combined to give the transonic potential (cid:13)ow equation
(2.5) (a2 u2)(cid:30) +(a2 v2)(cid:30) +(a2 w2)(cid:30)
xx yy zz
(cid:0) (cid:0) (cid:0)
2uv(cid:30) 2vw(cid:30) 2uw(cid:30) =0
xy yz xz
(cid:0) (cid:0) (cid:0)
in which u, v and w are the velocity components
(2.6) u=(cid:30) ; v =(cid:30) ; w =(cid:30)
x y z
The equation is elliptic at subsonic points for which q2 < a2, and hyperbolic at supersonic points for which
q2 >a2.
In a weak solution for an isentropic (cid:13)ow the appropriate jumps to be allowed conserve mass. Thus, at a
surface of discontinuity,
(2.7) n:((cid:26) q (cid:26) q )=0
1 1(cid:0) 2 2
where n is the normal to the surface, and the subscripts denote upstream and downstream conditions. At the
same time the introduction of a velocity potential ensures continuity of the tangential components of velocity.
Across such a jump the normal component of momentum is not conserved. The balance of momentum for the
complete(cid:13)owthenrequiresadragforceonthepro(cid:12)leequalandoppositetheforceonthejump. Thisrepresents
the wave drag according to the isentropic approximation (cf. [11]).
The boundary condition at the boundary is
@(cid:30)
(2.8) 0=q:(cid:23) =
@(cid:23)
where (cid:23) is the normal to the body surface. For lifting (cid:13)ows it is also necessary to impose the Kutta condition
which requiresthe circulationat eachsectionto be suchthat the velocityis (cid:12)nite at the sharptrailingedge. In
thethreedimensionalcasethisresultsinatrailingvortexsheetacrosswhichthereisajumpinpotential. Using
a simpli(cid:12)ed model, convection of the sheet will be ignored, and the sheet will be assumed to trail smoothly o(cid:11)
the trailing edge without rolling up. The conditions applied at the sheet are (cid:12)rst that the jump in potential
is constant along lines in the direction of the free stream, and second that the normal velocity is continuous
through the sheet. At in(cid:12)nity the velocity approaches a uniform far (cid:12)eld value except near the vortex sheet.
Thus (cid:30) , and it is necessary to work with a reduced potential from which the singularity at in(cid:12)nity has
! 1
been subtracted.
3
3 Coordinate Systems
It is di(cid:14)cult to satisfy the Neumann boundary condition (2.8) with su(cid:14)cient accuracy when the body surface
crosses the coordinate lines. This is particularly the case near the leading edge of the wing where the surface
has a high curvature and the (cid:13)ow is sensitive to small variations in the shape. The treatment is facilitated
by the use of a curvilinear coordinate system in which the body surface coincides with a coordinate surface.
In two-dimensional calculations such a system can be conveniently obtained by a conformal mapping. The
potential equation (2.5) becomes
(3.1) (a2 y2)(cid:30) 2uv(cid:30) +(a2 v2)(cid:30) +
xx xy yy
(cid:0) (cid:0) (cid:0)
(u2+v2)(uh +vh )=0
x y
where h is the modulus of the mapping function, x and y are new coordinates in the transformed plane, and u
and v are the velocity components in the new coordinate directions,
(cid:30) (cid:30)
x y
(3.2) u= ; v = :
h h
For two-dimensional (cid:13)ows which are subsonic at in(cid:12)nity, a most successful coordinate system, (cid:12)rst introduced
by Sells [12]forsubcritical calculations,is obtained bya mappingof the exteriorof the pro(cid:12)leontothe interior
of a circle. This coordinate system automatically bunches the mesh points in the sensitive regions near the
leading and trailing edges, and has been found to give high accuracy (cf. [5], [6]). Near the center of the circle
the modulus of the mapping equation approaches 1 where r is the radial coordinate, and the use of (cid:12)nite
r2
di(cid:11)erence formulas can lead to serious errors. It is convenient, therefore, to calculate the mapping onto the
exteriorofacircleandthentouseanexplicitinversion. Ifz isintheexteriorofthepro(cid:12)leandorintheexterior
of the circle, then the modulus of the mapping function may be expressed as
dz
h=
(cid:12)d(cid:27)(cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12)
where, if " is the included angle of the corner at the tra(cid:12)iling(cid:12) edge,
dz 1 1(cid:0)(cid:25)" N cn
(3.3) = 1 exp :
d(cid:27) (cid:18) (cid:0) (cid:27)(cid:19) (cid:26) (cid:27)n(cid:27)
nX=0
Themappingcoe(cid:14)cientsc canbecalculatedbyasimpleiterativeprocedure(cf. [13]). Ifthemappingfunction
n
is to be represented at 2K equally spaced grid points around the circle, it is convenient to take N = K so that
the real and imaginary parts of the power series reduce to trigonometric interpolation formulas for functions
speci(cid:12)ed at the grid points. This leads to favorable error estimates (cf. [14]). Evaluating the residue, the
condition that the mapping generates a closed pro(cid:12)le is seen to be
"
(3.4) c =1 :
1
(cid:0) (cid:25)
An airfoil with a wake of constant thickness is equally easily treated. It is simply necessaryto adjust the value
of c . If (cid:0) is the circulation, the far (cid:12)eld condition is (cf. [15])
1
cos((cid:18)+(cid:11)) (cid:0)
(3.5) (cid:30) + tan(cid:0)1[(1 M2)1=2tan((cid:18)+(cid:11))]
! r 2(cid:25) (cid:0)
where r and (cid:18) are polar coordinates for the interior of the circle and (cid:11) is the (cid:13)ow angle at in(cid:12)nity.
4
Ifthe(cid:13)owissupersonicatin(cid:12)nitythemappingtoacircleislesssatisfactorybecausewenowwishtodistinguish
between the conditions upstream at in(cid:12)nity where Cauchy data is to be provided, and downstream at in(cid:12)nity
where no condition should be imposed. A more convenient coordinate system is obtained by a mapping onto
the upper half-plane. Let w denote this plane. Then this mapping is easily derived from the mapping onto the
exterior of the circle by the additional transformation
1
(3.6) w=p(cid:27)+
p(cid:27)
The mapping from z to w, expressed with (cid:27) as a parameter, is then
N
dz dz d(cid:27) 2p(cid:27) c
n
(3.7) = = exp
dw d(cid:27)dw (1(cid:0)1=(cid:27))(cid:25)" (cid:26)nX=0(cid:27)n(cid:27)
For a cusped airfoil (" = 0) the singularities cancel at the trailing edge, leaving a smooth transformation.
Additional stretchings of the x and y coordinates can be used to map the half-plane to a rectangle. The (cid:13)ow
nowentersthe top boundary,splits, and leavesthroughthe twosides. Theairfoilcoversasegmentof the lower
boundary. Outside this segment corresponding points have to be identi(cid:12)ed as the two sides of a cut in the
z-plane, and the potential has to be continued acrossthe cut with a jump to accountfor the circulation. In the
far (cid:12)eld the mapping approaches a square root transformation
(3.8) z =w2
and a suitable reduced potential which is (cid:12)nite at in(cid:12)nity is obtained by setting
(3.9) G=(cid:30) (x2 y2)cos(cid:11) 2xysin(cid:11);
(cid:0) (cid:0) (cid:0)
wherexandyarecoordinatesinthetransformedplaneand(cid:11)isthe(cid:13)owangleatin(cid:12)nity. Ifthe(cid:13)owissubsonic
at in(cid:12)nity, the far (cid:12)eld condition becomes
(3.10) G=0 on the top boundary;
1
G= (cid:0) on the side boundaries;
(cid:6)2
where(cid:0)isthe circulation. Ifthe (cid:13)owissupersonicatin(cid:12)nity, itwillbeundisturbed upstreamof the bowwave,
so the condition is
@G
G=0; =0 on the top boundary .
@y
An economicalmethod of obtainingacoordinatesystem with similarcharacteristics,suggestedbyGarabedian,
issimplytoapplythesquareroottransformation(3.8)aboutapointjustinsidetheleadingedge. Thisgenerates
paraboliccoordinatesinthe physicalplane, whilethe pro(cid:12)leisunwrappedtoabump in thetransformedplane.
Unfortunately it is hard to satisfy the Neumann condition with su(cid:14)cient accuracy in the region of the leading
edgebecauseofthe rapidvariationofthe mappingmodulusnearthe singularityinsideit. Thisdi(cid:14)cultycanbe
overcomeby displacingthe y coordinatelines to follow the contourof the transformed pro(cid:12)le, so that it lies on
a coordinate line in a slightly nonorthogonal coordinate system (Figure 3.1). To avoid an unnecessary corner
in the coordinatescareshould be takento continuethe cut in the physicalplanesmoothly o(cid:11)the trailingedge.
It is impossible to removethe cornerdue to a (cid:12)nite wedgeangle at the trailing edge, but this is a feature to be
avoidedin any casebecauseit forcesthe appearanceof a stagnationpoint in the inviscid(cid:13)ow, leadingto larger
adversepressuregradients,andincreasingthe likelihoodof boundarylayerseparation. Themain disadvantages
of a sheared parabolic coordinate system of this type are the appearance of some extra terms in the equations
because the coordinates are not orthogonal, and less accurate resolution of the Kutta condition because there
is no concentration of mesh points near the trailing edge.
5
Figure 1:
In three-dimensional calculations the construction of a satisfactory coordinate system is one of the most
di(cid:14)cult parts of the problem. For studies of an isolated wing the use of sheared parabolic coordinates appears
to provide the best approach. For a tapered wing the use of a conformal transformation which varied in the
spanwisedirectionwouldaddnumerousextratermstotheequations,andtheresultingcoordinatesystemwould
in any case be nonorthogonal. Under a (cid:12)xed transformation independent of the spanwise coordinate, on the
other hand, the potential equation retains the relatively simple form
(3.11) (a2 u2)(cid:30) +(a2 v2)(cid:30) +h2(a2 w2)(cid:30)
xx yy zz
(cid:0) (cid:0) (cid:0)
2uv(cid:30) 2hvw(cid:30) 2huw(cid:30) +(u2+v2)(uh +vh )=0
xy yz xz x y
(cid:0) (cid:0) (cid:0)
Thus if the wing hasastraightleadingedge, which is the mostusual case,the squareroottransformation(3.8)
can conveniently be used to unwrap the sections about a singular line just inside the leading edge, so that the
wing surface becomes a shallow bump above a coordinate plane formed by unfolding the two sides of a cut
behind the singular line. Then, as in the two-dimensional case, the bump can he removed by displacing the
coordinatesurfacessothat they become parallelwith the wing surface. The trailingvortexsheet is assumedto
lie along the cut, so that it is also split by the transformation. A complication is caused by the continuation of
the cut beyond the wing tips. Points on the two sides of the cut must be identi(cid:12)ed as representing the same
physicalpoint. Alsoaspecialformoftheequationsmustbeusedatpointsonthe‘fold’alongthecontinuationof
the singular line, where h vanishes. At these points the equation to be satis(cid:12)ed reduces to the two-dimensional
Laplace equation in the x and y coordinates.
Analyticalexpressionsforthefar(cid:12)eldpotentialhavebeengivenbyKlunker[16]. If,however,the coordinates
are stretched to in(cid:12)nity, the use of these can be avoided because at in(cid:12)nity the square root transformation
collapses the region in(cid:13)uenced by the vortex sheet to a single line of mesh points containing the sheet.
Thecaseofayawedwing,asproposedbyR.T.Jones,istreatedbykeepingthecoordinatesystem(cid:12)xedtothe
wing, with the x,y-planes normal to the leading edge, and rotatingthe (cid:13)ow at in(cid:12)nity through the appropriate
yaw angle, as illustrated in Figure 3.2. The vortex sheet then has to be tracked in the stream direction behind
the trailing edge and downstream tip.
4 Simple and Rotated Di(cid:11)erence Schemes
Thenumericalformulationoftheproblemrequirestheconstructionofanappropriatesetofdi(cid:11)erenceequations.
Thesearethentobesolvedbyiteration. FollowingtheideaadvancedbyMurmanandCole[3],theplanofattack
is to distinguish between the regions of subsonic (cid:13)ow, where the governing equation is elliptic, and supersonic
6
Figure 2: Yawed wing.
(cid:13)ow,whereitishyperbolic,andtouseadi(cid:11)erenceschemewhichissensitivetothe localtype. Forthis purpose
the velocities are (cid:12)rst determined, using central di(cid:11)erence formulas to approximate the (cid:12)rst derivatives of the
potential. Then at elliptic points central di(cid:11)erence formulas are used for all second derivatives. At hyperbolic
pointsstreamwisesecondderivativesareapproximatedbydi(cid:11)erenceformulasretardedintheupstreamdirection.
This yields an implicit system of equations with a region of dependence which correctly re(cid:13)ects the region of
dependence of the (cid:13)ow, and leads to a truncation error which has the e(cid:11)ect of an arti(cid:12)cial viscosity. It is
important to use an implicit di(cid:11)erence scheme, because near the sonic line the characteristics become almost
perpendicular to the velocity, with the result that the mesh spacing of an explicit scheme would have to be
reducedtozerointhedirectionof(cid:13)owtoobtainthecorrectregionofdependence. Thearti(cid:12)cialviscosityensures
thedesiredentropyinequality,excludingexpansionshocksbutpermittingtheappearanceofcompressionlayers.
If the (cid:13)ow is aligned with one of the coordinates, it is relatively easy to implement these ideas. Considering
the two-dimensionalcase(3.1) for convenience,let x be the coordinatein the (cid:13)owdirection. Then if the (cid:13)owis
supersonic at the mesh point [i(cid:1)x;j(cid:1)y], we take
(cid:30)ij 2(cid:30)i(cid:0)1;j +(cid:30)i(cid:0)2;j
(4.1) (cid:30) = (cid:0) ;
xx (cid:1)x2
(cid:30)i;j+1 (cid:30)i;j(cid:0)1 (cid:30)i(cid:0)1;j+1+(cid:30)i(cid:0)1;j(cid:0)1
(cid:30) = (cid:0) (cid:0)
xy
2(cid:1)x(cid:1)y
(cid:30)i;j+1 2(cid:30)ij +(cid:30)i;j(cid:0)1
(cid:30) = (cid:0)
yy (cid:1)y2
Thus we obtain a simple implicit scheme in the hyperbolic zone which can be solved line by line. The
schemeissecondorderaccuratein yand (cid:12)rstorderaccuratein x. Thedominanttruncationerrorsarisingfrom
the terms in (cid:30) and (cid:30) are (cid:1)x(u2 a2)(cid:30) and (cid:1)xuv(cid:30) . These represent a positive arti(cid:12)cial viscosity
xx yy xxx xy
(cid:0)
7
provided that u2 > a2. When the (cid:13)ow is not perfectly aligned with the coordinate system, there exist points
in the hyperbolic region at which u2 <a2 <u2+v2. These are a source of di(cid:14)culty. On the one hand, a von
Neumanntestindicates,andexperiencecon(cid:12)rms,thatcentraldi(cid:11)erenceformulasareunstableinthehyperbolic
region. Ontheotherhand, theycoordinatelinenowliesbehindoneofthecharacteristicsfromthecenterpoint
[i(cid:1)x;j(cid:1)y], so that use of the retarded di(cid:11)erence formulas leads to a scheme which does not have the correct
region of dependence. This is re(cid:13)ected in the fact that the arti(cid:12)cial viscosity (cid:1)x(u2a2)(cid:30) is now negative,
xxx
and the existence of a marching instability in the x direction is con(cid:12)rmed by a von Neumann test.
Nevertheless, the scheme has proved extremely satisfactory in practice, provided that the supersonic (cid:13)ow is
con(cid:12)ned to a zone near the body in which the velocity is more or less aligned with the curvilinear coordinate
system. The presence of a narrow band of points with negative numerical viscosity near the sonic line is not
enough to upset the calculation. Instability can be anticipatedas the mesh spacingis shrunk to zero,but more
than adequate accuracy can be obtained well before its onset. Following Garabedian and Korn [5], one may
also improve the accuracy by partially cancelling the truncation errorsthrough the addition of (cid:12)nite di(cid:11)erence
formulas representing (cid:1)x(cid:30) and (cid:1)x(cid:30)xxy .
xxx
When the (cid:13)owissupersonicatin(cid:12)nity, andalsointhree-dimensionalcalculations,asimpledi(cid:11)erencescheme
of this type is not so satisfactory. Curvilinear coordinate systems which allow an accurate treatment of the
boundary condition are generally not well aligned with the (cid:13)ow. We need, therefore, a coordinate invariant
di(cid:11)erence scheme. Considering (cid:12)rst the two-dimensional case, the principal part of (3.1) can be written in the
canonical form
(4.2) (a2 q2)(cid:30) +a2(cid:30) =0;
ss nn
(cid:0)
where s and n are coordinates in the local stream and normal directions. Since the direction cosines are u=q
and v=q, (cid:30) and (cid:30) can be expressed locally in terms of the actual coordinates as
ss nn
1
(4.3) (cid:30) = (u2(cid:30) +2uv(cid:30) +v2(cid:30) )
ss q2 xx xy yy
1
(cid:30) = (v2(cid:30) 2uv(cid:30) +u2(cid:30) )
nn q2 xx(cid:0) xy yy
It is now easy to devise a proper di(cid:11)erence scheme for points in the supersonic zone. By using di(cid:11)erence
formulasretardedinboththexandydirectionsforallcontributionsto(cid:30) ,andcentraldi(cid:11)erenceformulasfor
ss
all contributions to (cid:30) , we obtain a ‘rotated’ di(cid:11)erence scheme which reduces to the Murman scheme when
nn
the (cid:13)ow is aligned with either coordinate direction. At points in the subsonic zone central di(cid:11)erence formulas
are used in the conventional manner.
A positive arti(cid:12)cialviscositywith a magnitude proportional to q2 a2 is introduced at all supersonicpoints.
(cid:0)
This enforces the correct directional property on the solution. At the same time the presence of points ahead
of the y coordinateline in the formulas for (cid:30) , ensures that the region of dependence of the di(cid:11)erence scheme
nn
alwayscontainsthe regionof dependenceof the di(cid:11)erentialequation. Inthe subsoniczonethe schemeissecond
order accurate. In the supersonic zone it is second order accurate in the term representing (cid:30) , but only (cid:12)rst
nn
order accurate in the term representing (a2 q2)(cid:30) , forcing the use of rather (cid:12)ne mesh to obtain an accurate
ss
(cid:0)
solution.
5 Iterative Solution: Analysis Using Arti(cid:12)cial Time
While the coordinateinvariant form of the rotated di(cid:11)erence scheme assuressatisfaction of the correctentropy
inequality,italsoincreasesthedi(cid:14)cultyofactuallysolvingthedi(cid:11)erenceequations. Theuseofcentraldi(cid:11)erence
formulas in both coordinate directions for the contributions to (cid:30) prevents the use of a simple marching
nn
procedure to solve the equations in the hyperbolic region a line at a time, as in the Murman scheme. Instead,
8
if one advances in the x direction, the correction at the point [i(cid:1)x;j(cid:1)y] has to be calculated using ‘old’ values
for the forward points [(i+1)(cid:1)x;(j +1)(cid:1)y];[(i+1)(cid:1)x;j(cid:1)y] and [(i+1)(cid:1)x;(j 1)(cid:1)y]. Thus the equations
(cid:0)
in the hyperbolic region have to be solved by iteration, just as in the elliptic part of the (cid:13)ow.
Such a process has to be carefully controlled to ensure stability and convergence. It is helpful to regard the
iterations as successive levels in an arti(cid:12)cial time coordinate. This leads to a theory for equations of mixed
type in the spirit of Garabedian’s theory for elliptic equations in [17]. The essential idea is that the iterative
procedure, regarded as a (cid:12)nite di(cid:11)erence scheme for a time dependent process, should be consistent with a
properly posed initial value problem which is ‘compatible’ with the steady state equation; that is to say, that
the solution of the initial value problem should reach an equilibrium point which represents a solution of the
steady state equation. Then, if the di(cid:11)erence scheme is also stable, it should convergeto the desired solution.
Let (cid:1)t denote the time step, and let updated values of (cid:30) at any circle of the calculation be denoted by the
superscript +. Then the typical central di(cid:11)erence formula for a second derivative is
(cid:30)+ (1+r(cid:1)x)(cid:30)+ (1 r(cid:1)x)(cid:30) +(cid:30)
(cid:30) = i(cid:0)1;j (cid:0) ij (cid:0) (cid:0) ij i+1;j
xx (cid:1)x2
where the updated value is used on one side, the old value on the other side because the updated value is not
yet available, and a linear combination of the two is used at the center point. In the time dependent system
this formula may be interpreted as representing
(cid:1)t
(cid:30) ((cid:30) +r(cid:30) )
xx xt t
(cid:0) (cid:1)x
We must consider,therefore,a time dependent equationwhich containsmixed space-timederivatives. Dividing
through by a2, its principal part will have the form
(5.1) (M2 1)(cid:30) +2(cid:11)(cid:30) (cid:30) +2(cid:12)(cid:30) =0
ss st nn nt
(cid:0) (cid:0)
where M is the local Mach number q=a, and the values of (cid:11) and (cid:12) depend on the combination of new and old
values actually used in the di(cid:11)erence formulas. Introducing a new time coordinate,
(cid:11)s
(5.2) T =t +(cid:12)n
(cid:0) M2 1
(cid:0)
(5.1) becomes
(cid:11)2
(5.3) (M2 1)(cid:30) (cid:30) (cid:12)2 (cid:30) =0
(cid:0) ss(cid:0) nn(cid:0)(cid:18)M2 1 (cid:0) (cid:19) TT
(cid:0)
If the (cid:13)ow is locally supersonic, it can be seen that T is spacelike and either s or n is timelike. Since s is the
timelike direction in the steady state problem, we require
(5.4) (cid:11)2 >(cid:12)2(M2 1)
(cid:0)
to ensure that it is also the time-like direction in the unsteady problem. This indicates the need to add terms
representing mixed space-time derivatives to the retarded di(cid:11)erence formulas to compensate for their presence
in the central di(cid:11)erence formulas.
The characteristicsof the unsteady equation (5.1) satisfy
(5.5) (M2 1)(t2+2(cid:12)nt) 2(cid:11)st+((cid:12)s (cid:11)n)2 =0
(cid:0) (cid:0) (cid:0)
Provided that (cid:11) is positive, the cone of dependence lies on the upstream side of the n,t-plane and behind the
s,n-plane, touching along the line
(cid:12)
n= s
(cid:11)
9
Figure 3: Characteristic cone of equivalent time dependent equation.
as illustrated in Figure 5.1. Thus the region of dependence at the current time level is reduced to a single
line, and an iterative scheme with the correct region of dependence can be obtained simply by ensuring that
the points surrounding this line have already been updated. If (cid:11) is negative, the cone of dependence is in
the downstream direction. Thus to ensure the proper region of dependence it is necessary to control both the
magnitude and the sign of the coe(cid:14)cient of (cid:30)
st
A von Neumann test of local stability leads to an additional guideline. Let the di(cid:11)erential equation be
approximated at the mesh point [i(cid:1)x;j(cid:1)y] by the general di(cid:11)erence equation
(a (cid:30) b (cid:30)+ )=0
pq i+p;j+q (cid:0) pq i+p;j+q
Xp;q
Substituting the formula
(cid:30)=Gkeimxeiny
for (cid:30) at the k-th time interval, the growth factor G is given by
a ei(p(cid:24)+q(cid:17))
G= p;q pq
P p;qbpqei(p(cid:24)q(cid:17))
P
where
(cid:24) =m(cid:1)x; (cid:17) =n(cid:1)y:
For small values of (cid:24) and (cid:17), G may be expanded as
A +A i(cid:24)+A i(cid:17)+A (cid:24)2+A (cid:24)(cid:17)+A (cid:17)2+:::
00 10 01 20 11 02
G=
B +B i(cid:24)+B i(cid:17)+B (cid:24)2+B (cid:24)(cid:17)+B (cid:17)2+:::
00 10 01 20 11 02
where consistency with the steady state di(cid:11)erential equation (3.1) requires that
A =B ; A =B ; A =B
00 00 10 10 01 01;
a2 u2 2uv a2 v2
B A = (cid:0) ; B A = (cid:0) ; B A = (cid:0)
20(cid:0) 20 (cid:1)x2 11(cid:0) 11 (cid:1)x(cid:1)y 02(cid:0) 02 (cid:1)y2
Since (cid:24) and (cid:17) can be chosen so that
(cid:24)2 2uv(cid:24)(cid:17) (cid:17)2
(a2 u2) 2 +(a2 v2)
(cid:0) (cid:1)x2 (cid:0) (cid:1)x(cid:1)y (cid:0) (cid:1)y2
10