1. science
2. mathematics
3. geometry

AI AA-85-0435
Finite Volume Solution of the
Two-Dimensiona I Euler Equations
on a Regular Triangular Mesh
A.Jameson and D. Mavriplis
Princeton University, Princeton, N.J.
AlAA 23rd Aerospace Sciences Meeting
January 14-17,1985/Reno, Nevada
uc"
For permission to copy or republish, Contact the American Institute of Aeronautm and Astronautics
1633 Broadway, New York. NY 10019
FINITE VOLUME SOLUTIOK OF THE TWO-DIMENSIONAL
EULER EQUATIOXS OX A REGULAR TRIANGULAR MESH
by
A. Jameson and D. Mavriplis
DeDartment of Mechanical and Aerospace Engineering
Princeton University
Princeton, New Jersey 08544
Abstract
engine mcelle/pylon combinations to approximate
the entire aircraft. The major obstacle in
solving the flow problem over these increasingly
complex geometries is the generation of smoothly
varying meshes about the body configuration with
appropriate cell distributions to calculate the
flow variables over the entire body surface with
sufficient accuracy. Various approaches to overcome this obstacle have been proposed, including
local mesh enrichment',
and hybrid meshes5, for
example a global +mesh with imbedded 0-meshes
around the components of interest. Another
approach is to use alternative cell geometries
such as triangles or polygons. The flow solution
about an entire aircraft will probably require the
use of all these approaches simulceneously. This
work is concerned with the use of triangdar elements, which allow an added degree of flexibility
in fitting complex geometries. Recently ~ e l z 6
obtained the transonic potential flow solution on
regular and irregular triangular grids in two
dimensions, and showed how these could be applied
to complex geometries not readily amenable to
quadrilateral meshes.
The two dimensional Euler equations have been
solved on a regular triangular grid using the
finite volume approach. By careful construction
of the dissipative terms, the scheme is designed
to be second order accurate in space, provided the
grid is smooth, except in the vicinity of shocks,
where it behaves as first order accurate. The
multigrid technique has been applled to triangular
meshes to enhance the convergence of the scheme.
In its present form, the accuracy and convergence
rate of the triangle code are comparable to that
of the quadrilateral mesh code of 3ameson.l
1. Introduction
w
In the last decade, transonic flow calculations using the potential flow assumption have
become fairly widespread and have proved accurate
and cheap enough to be regularly used as a design
tool. However, in the transonic regime the potential flow assumption is not strictly correct for
it neglects the entropy and vorticity production
associated with shock waves.
In this work, the Euler equations have been
solved on a regular triangular grid about a simple
airfoil geometry in two dimensions. The solution
scheme is based upon that of Jameson, Schmidt and
Turke12 and convergence has been accelerated
by several means including the use of a multigrid
scheme. It is concluded that the triangular mesh
scheme can be made competitive with the quadrilateral mesh scheme. Thus, triangular grids
should prove useful in the application to complex
geometries.
To describe inviscid transonic flow correctly,
the Euler equations rmst be solved. The numerical
s01ution of the Euler equations should provide an
accurate prediction of the location and strength
of a shock and the associated wave drag.
Furthermore, the solution should be capable of
predicting rotational flow problems such as the
slip stream of a prop-fan flowing over a wing, or
the flow over aircraft components in the vortex
wake downstream of a wing.
2. Mesh Generation
b#
Recent developments have substantially reduced
the cost of Euler calculations. The finite volume
formulation of Jameson, Schmidt and Turke12
which uses a Runge-Kutta type time stepping scheme
has proved to be robust and accurate. The rate of
convergence has also been greatly enhanced by
various means, including the use of a rmltigrid
scheme. It has in fact been shown in a recent
paper3 that this scheme is now sophisticated
enough to be applied to complete
winglhodyltail-planelfin combinations. The next
step is to apply the scheme to geometries of
increased complexity, such as the addition of
-1
A two-dimensional regular triangular mesh is
generated about an airfoil in similar fashion to
the generation of the cormnon quadrilateral C-mesh.
The airfoil is first conformally mapped to a near
circle by a Karma" Trefftz transformation. A
family of similar near-circles with increasing
radii can now be generated, together with radial
coordinate lines at regular azimuthal intervals.
In this form, the radial and nearcircle coordinate line intersections form the nodes of a
quadrilateral grid. Each quadrilateral may be
split into two triangles by drawing the diagonal
@ A m . W R l t Of A m ~ . n t l n .Dd
-mtks.Ir..5?85.Mritbllrrmcd.
1
j o i n i n g t h e lower l e f t hand c o r n e r t o t h e upper
r i g h t hand c o r n e r of t h a t cell. I n t h i s manner an
o r t h o g o n a l q u a d r i l a t e r a l mesh is transformed i n t o
a mesh of r i g h t a n g l e t r i a n g l e s .
I f t h e mesh
nodes on e v e r y second n e a r c i r c l e are d i s p l a c e d by
h a l f a mesh w i d t h , a symmetric t r i a n g u l a r mesh is
obtained.
The o r i g i n a l mapping is f i n a l l y
i n v e r t e d t o produce a r e g u l a r t r i a n g u l a r mesh
about t h e a i r f o i l .
d i s c r e t i z a t i o n of t h e E u l e r e q u a t i o n s l e a d s t o a
s y s t e m of coupled o r d i n a r y d i f f e r e n t i a l e q u a t i o n s
+ Q(vi)
%sini)
On a t r i a n g u l a r mesh, t h e r e l a t i o n between t h e
4.
The v a r i a b l e s t o be determined are t h e
p r e s s u r e , d e n s i t y , C a r t e s i a n v e l o c i t y components,
t o t a l energy and t o t a l e n t h a l p y , denoted by p, P ,
u , v , E and H r e s p e c t i v e l y . Since f o r a p e r f e c t
g a s we have t h e r e l a t i o n s
P
+'
2
dik =
II
n
w dxdy +
1 (fdy
afi
- gdx)
= 0
TE ( 4 ) [-Asit i + -Sk
1lV
Atk
an..
-
2
V wkl
1 Wk -
3Wi
.
and is of o r d e r (Ax)
Thus t h e d i s s i p a t i v e
o p e r a t o r D is of o r d e r (Ax)3 r e l a t i v e t o t h e conv e c t i v e o p e r a t o r Q, and the accuracy of t h e scheme
is p r e s e r v e d .
In r e g i o n s of s t r o n g p r e s s u r e grad i e n t s s u c h as near a shock wave. t h i s d i s s i p a t i o n
i s i n s u f f i c i e n t t o s u p p r e s s o s c i l l a t i o n s and we
must add l a r g e r term. The d i k above are t h u s t o
be augmented by
PUV
PVH
PUH
I
k=l
x and Y
2
PE
w.
3
(1)
V2Wi =
where n is a f i x e d area With boundary.
are C a r t e s i a n c o o r d i n a t e s and
(3)
The S l a t term is p r o p o r t i o n a l t o the size of the
c e l l f a c e i k and scales a p p r o p r i a t e i y w i t h t h e
time d e r i v a t i v e in e q u a t i o n ( 3 ) . V is t h e undiv i d e d 4 p o i n t L a p l a e i a n o p e r a t o r on a twod i m e n s i o n a l t r i a n g u l a r mesh:
where I is t h e r a t i o of s p e c i f i c h e a t s , we need
o n l y solve f o r t h e f o u r v a r i a b l e s P , P U , PV, and
pE. These values are determined by s o l v i n g t h e
E u l e r e q u a t i o n s , which i n i n t e g r a l form read:
a
= 0
where d i k is t h e f l u x across t h e f a c e d e l i m i t i n g
c e l l i and i t s neighbour k.
H = E + E
P
(u2+v2),
-
where D(w )must have a s u i t a b l e form t o damp o u t
I n keeping w i t h t h e
a l l undesired o s c i l l a t i o n s .
f i n i t e volume f o r m u l a t i o n , D(w ) is c o n s t r u c t e d as
a summation of f l u x e s across t k e cell f a c e s :
The s o l u t i o n procedure f o r t h e E u l e r e q u a t i o n s
on a t r i a n g u l a r mesh c l o s e l y f o l l o w s t h a t proposed
by Jameson, Schmidt and Turke12 f o r q u a d r i l a t e r a l
meshes.
-
Dissipation
+ Q(wi)
&tSiwi)
D i s c r e t i z a t i o n of t h e E u l e r E q u a t i o n s
(Y-l)P
E q u a t i o n (1) can be a p p l i e d t o each i n d i v i d u a l
t r i a n g u l a r c e l l . For t h e rmomentum, f o r example,
it yields
3
a
(SiPul)
+
(Qkpuk + byk P,) = 0
(2)
ik
si
[2
Ati
E
1
+ k1-A' t k [wi -
w I
k
k=l
where
where S is t h e c e l l a r e a , Qk is t h e f l u x v e l o c i t y
a c r o s s s i d e k of t h e cell,
Q, = byk
\ -
A\
w
As is t h e case w i t h q u a d r i l a t e r a l meshes, c e l l
c e n t e r e d E u l e r schemes on t r i a n g u l a r meshes a l l o w
d e c o u p l i n g between odd and even c e l l s . In o r d e r
t o p r e v e n t t h i s d e c o u p l i n g e x t r a d i s s i p a t i v e terms
are needed. Thus e q u a t i o n ( 2 ) is r e p l a c e d by
depending on an add or even row and column number,
from which t h e l o g i c f o r t h e numbering can be
deduced. The numbering scheme is d e p i c t e d i n
F i g u r e 1.
E
(2)
where Q is t h e d i s c r e t e approximation t o t h e flux
i n t e g r a l i n (1).
c e l l numbering and g r i d p o i n t numbering can be
complicated, 8s may be t h e r e l a t i o n between a part i c u l a r c e l l number and t h a t of i t s n e i g h b o r s .
For t h e s p e c i a l case of a r e g u l a r t r i a n g u l a r mesh,
four d i f f e r e n t t y p e s of c e l l s can be i d e n t i f i e d .
3.
f a r each c e l l i
= 0
s h o u l d be of o r d e r 1 near a shock and
of o r d e r (Ax)2 in regions of smooth flow where t h e
2nd o r d e r accuracy of t h e scheme mst be preserved.
One method of a c h i e v i n g t h i s i s t o t a k e
vk
ci:)
p r o p o r t i o n a l t o t h e undivided 4 p o i n t
2
L a p l a c i a n of t h e p r e s s u r e , V p i . This, however,
was not found t o p r o v i d e s a t i s f a c t o r y b e h a v i o r .
I f t h e boundary is a l i g n e d v i t h one of t h e coord i n a t e lines, p r e s s u r e is e x e r t e d on i t only by
e v e r y 2nd , c e l l (see Fig. 2). ,Thus we d e f i n e v i as
and A$,
by are t h e i n c r e m e n t s of x and y a l o n g
s i d e k of tke cell. The flow v a r i a b l e s such as
P U ~a r e t a k e n as t h e a v e r a g e of t h e two cell
v a l u e s on e i t h e r s i d e of t h a t f a c e . On a smooth
g r i d , t h i s r e d u c e s t o c e n t r a l d i f f e r e n c i n g and i s
second o r d e r a c c u r a t e in space. I n t h i s form, t h e
2
upi
where wn and wn+l are t h e v a l u e s a t t h e beginning
and end of t h e nth t i m e s t e p . The s t a n d a r d v a l u e s
f o r the c o e f f i c i e n t s are
W
where t h e s u b s c r i p t s j have been o m i t t e d s i n c e
t h e y are c o n s t a n t . In t h i s c a s e , j = c o n s t a n t
d e f i n e s t h e c o o r d i n a t e l i n e of i n t e r e s t . Now Uik
i s d e f i n e d as t h e maximum v a l u e of v i f o r t h e f o u r
cells on e i t h e r s i d e of f a c e i k s l a n g t h e coordinate line j
constant (i.e. i f face i k delimits
c e l l s i,j and i + l , j , v i k is t h e uximum of t h e s e t
vi-3
t h r o u g h vi+&).
For each f a c e i k , one of t h e
3 c o o r d i n a t e d i r e c t i o n s is a l i g n e d w i t h t h i s f a c e ,
and t h u s o n l y 2 d i r e c t i o n s c o n t a i n components normal t o t h e f a c e i k . Hence we d e f i n e
T h i s p a r t i c u l a r scheme is used because i t is
well s u i t e d t o drive t h e r m l t i g r i d a l g o r i t h m .
To
a c c e l e r a t e convergence, A t is t a k e n as t h e maximum
p e r m i s s i b l e l o c a l time s t e p s e p a r a t e l y i n e a c h
c e l l . This g u a r a n t e e s that d i s t u r b a n c e s will be
e x p e l l e d t o t h e Outer boundary i n a f i x e d number
of s t e p s p r o p o r t i o n a l t o t h e number of cells between the i n n e r and o u t e r boundary.
-
(2)
cik
-
6.
9, = 1,2,3
8
k2 =x(uik)
The boundary c o n d i t i o n s are t h o s e d e s c r i b e d by
Jameson, Schmidt and T u r k e l l and Jameson and
Baker3, and w i l l o n l y be o u t l i n e d b r i e f l y .
where lreik is t h e u i k c a l c u l a t e d along t h e m o r d i n a t e line e and k2 is a c o n s t a n t c o e f f i c i e n t .
A t t h e i n n e r boundary t h e r e is no f l u x through
t h e w a l l . However, t h e v a l u e s of t h e p r e s s u r e at
t h e w a l l are needed t o c a l c u l a t e t h e c o n t r i b u t i o n s
bXk pk and Ayk pk t o t h e mOmentUm e q u a t i o n .
The
wall p r e s s u r e can be e x t r a p o l a t e d from t h e
p r e s s u r e a t t h e c e n t e r of t h e boundary cell by use
of t h e e q u a t i o n s g i v e n in r e f e r e n c e 2 :
In t h e v i c i n i t y of a shock wave, t h e 4 t h d i f f e r e n c e s tend t o produce o v e r s h o o t s , hence t h e y
a r e t u r n e d o f f by d e f i n i n g
E(4)
ik
Boundary C o n d i t i o n s
w i t h P = k omitted
= max(0.
k4
-
E
(2))
ik
where k4 is a c o n s t a n t c o e f f i c i e n t .
The f i n a l form of t h e d i s s i p a t i v e f l u x e s i s t h u s :
d i k = - [ -1+ - si
)
2 Ati
k
'
Atk
5.
where X and Y r e p r e s e n t c o o r d i n a t e l i n e s a l i g n e d
and normal t o t h e boundary and x and y are t h e
r e g u l a r C a r t e s i a n c o o r d i n a t e s . This e q u a t i o n
r e l a t e s t h e normal p r e s s u r e g r a d i e n t a t t h e wall
due t o c e n t r i f u g a l f o r c e s t o t h e c u r v a t u r e of t h e
boundary.
-
A t t h e e x t e r i o r boundary we wish t o minimize
t h e r e f l e c t i o n of o u t g o i n g d i s t u r b a n c e s . C o n s i d e r
t h e flow normal t o t h i s boundary.
Assuming i t t o
be l o c a l l y one-dimensional, we i n t r o d u c e t h e f i x e d
and e x t r a p o l a t e d Riemann i n v a r i a n t s :
Time S t e p p i n g
A f i v e s t a g e time s t e p p i n g scheme is used
where', f o r economy t h e d i s s i p a t i v e o p e r a t o r is
e v a l u a t e d only i n t h e f i r s t 2 s t a g e s .
Since t h e
cell volume Si is c o n s c a n t . e q u a t i o n ( 3 ) may be
w r i t t e n as
dwi 1
dt + Si (Q(w,) - D(wi)l
= 0
c o r r e s p o n d i n g t o incoming and o u t g o i n g charact e r i s t i c s . The normal v e l o c i t y and l o c a l speed of
sound m y t h u s be determined by
We t h e n advance i n time by t h e scheme
( 0 ) = wn
W
1
n --(Re
2
-
q
c
-%Re
+
R,)
- R,)
Two o t h e r independent c o n d i t i o n s a r e needed t o
complete t h e d e f i n i t i o n of the o u t e r b u n d a r y cond i t i o n . These are g i v e n by t h e v a l v e s of tangent i a l v e l o c i t y and e n t r o p y . For an o u t f l o w
boundary r h e s e are e x t r a p o l a t e d from t h e i n t e r i o r
v a l u e s whereas f o r an i n f l o w boundary t h e y are set
equal t o t h e i r f r e e s t r e a m v a l u e s .
3
.*
performing several Jaeobi iterations on this
system. In practice E is assigned the value 112
and only 2 iterations are required to establish
the approximate values of the smoothed residuals
necessary for accelerating convergence of the time
stepping scheme.
7. Convergence Acceleration
Various devices have been employed to accelerate the convergence of the Euler equations to a
steady state. One is the use of a local time step
calculated from local flow properties, such that
the scheme operates everywhere at its stability
limit. The Euler equations themselves may also be
modified far increased convergence, provided the
steady state solution is not altered. Since, in
the steady state H = tl, throughout the flowfield,
it may be argued that we need not integrate the
energy equation, but simply keep the enthalpy
constant at its free-stream value. The pressure
may then be calculated from the values of density
and velocity. Gm the other hand, we may integrate
the energy equation and use the difference H-H,
as a forcing function to accelerate convergence.
apu(H-H,),
Then artificial terms ap(H-H,),
apv(H-Jb,) and a p ( H - L ) are added to the mass,
momentum and energy equations respectively, where
o is B constant. Enthalpy damping of this type
has proved effective i n practice2.
irvi
ii) Kultigrid Scheme
The idea of the multigrid scheme is to use
corrections calculated on successively coarser
grids to improve the solution on a fine grid. The
motivation stems from the observation that the
time stepping scheme is capable of rapidly d i m nating high frequency errors, but is very slow at
reducing lower frequency ot global errors. Since
the highest frequency errors which are visible to
a particular grid are of the order of one mesh
width, time stepping on progressively coarser
grids will have the effect of rapidly elirmnating
a distinct bandwidth of ertors on each grid.
The 5 stage scheme is used to drive the
multigrid algorithm because it has excellent high
frequency damping properties. For a regular
triangular mesh, coarser grids may be constructed
by doubling the mesh width and assembling groups
of 4 small triangles to form one large triangular
cell (Fig. (3)).
The flow variables are transferred to the coarser grid by the rule
i) Residual Averaging
The maximum time step which may be taken is
limited by the Courant-Fredrichs-Lewy condition
which states that the domain of dependence of the
discretized equations m s t at least contain that
of the original differential equation. To ease
this testriction we may perform averaging on the
residuals to increase the support of the scheme.
If the residual at cell i is
we replace it by
(O)
W2h
-1
(Sh Wh)/SZh
where the sabscripts denote the mesh spacing. S is
the cell area, and the sum i s over the 4 triangles
composing each cell on the coarse grid. This =le
conserves mass, mmentum and energy. We may similarly sum the fine grid residuals to obtain
R Zh
( O ) = 1 I$, (wh)
q:
Our objective i s to find a correction 6w such
that
where E is a constant and V 2 is the 4 point undivided Laplaeian operator. With this procedure it
vas found that the Courant number At/Ax could be
increased by about one half, but little improvement in convergence was observed. It has been
shown that on quadrilateral meshes performing the
smoothing implicitly leads to infinite support,
and that the Courant-Fredrichs-Levy condition no
longer imposes any limit on the Courant number.
Following this approach, we set
1 Rh(w
+ 6w) = 0
Remembering that the residual is a nonlinear function of w, we approximate this as
RZh(wZh + 6w) RZh(wzh) + RZh
(0) = 0
(4)
-
a
)R
::
where (
7
6w has)
been replaced by the
corresponding difference an the coarse grid.
Since the left hand side of equation ( 4 ) is the
quantity we wish to minimize, we advance in time
on the coarse grid a8 follows
where the barred values indicate new smoothed
residuals. This may be rewritten as
for the qth stage of the time stepping scheme.
Thus, time stepping on the coarse grid proceeds 8 s
on the fine grid except for the addition of the
forcing function
The coefficient matrix of this system of equations
contains 4 no"-zero elements in each row, and eannot be exactly inverted by an inexpensive
algorithm. However, for values of (. in the range
of interest ( i . e . , usually less than 1) this
matrix is strongly diagonallxdominant. Thus an
approximate solution to the Ri may be obtained by
4
.*rar
Z?
e
attack and Mach number of .75. Horawetz. ;heorem7
guarantees that any shock-free solution in transonic flow is an isolated point. The KORN airfoil
is designed by the hodograph method to attain such
a shock-free solution at the above conditions.
Thus the shock-free solution obtained numerically
and shown i n Figure 11 is a strong indication of
the accuracy of the scheme. Finally, a subcritical lifting case is shown in Figure 12 where the
drag 16 seen to be virtually zero, as expected.
The lift and drag coefficients for all these cases
agree to within 3 decimal places with those caleulated using FL052.
This may be repeated on successively coarser
grids. Om returning from the coarse to the fine
grid, if is found that each fine grid cell center
lies on a straight line which joins two coarse
grid cell centers, and the corrections are therefore linearly interpolated between the two coarse
grid cell centers. This is illustrated in Figure
3.
An effective mltigrid strategy is to calculate 1 rime step on each grid when passing to
coarser grids, and simply to interpolate the flow
variables without time stepping when returning to
finer grids.
FOK subcritical eases, where no shock
waves are present, the scheme should be second
order accurate provided that the grid is smooth
enough. Furthermore, Since the flow will be irrotational, the drag should be zero. and the calculated value of CD m y provide an estimate of the
discretization error. The subcritical lifting
case of a NACA 0012 airfoil at Mach number 0.5 and
3O angle of attack has been run for variom grid
sizes and the value6 of CL and CD are given in
Table 1 as a function of grid size. These results
verify the second order accuracy of the scheme.
In practice, computations begin by imposing a
uniform flow field over the entire grid. At time
t = 0 the airfoil is impulsively introduced by
suddenly requiring that the boundary condition at
its surface be satisfied. The Euler equations are
then integrated in time until a steady state is
reached. Convergence can thus be monitored by the
average size of the residuals in the flow field.
At each stage in a single time step the convective
operator Q is recalculated and the dissipative
operator D is either recalculated or taken as that
of the previous stage. The resulting residuals
are then smoothed and the flow variables are
updated by the time stepping scheme. After all
stages of the time step are completed. the flow
variables are modified by introducing the enthalpy
damping terms.
w
The convergence Kate of the triangle code for
supercritical case on a coarser 64x16 mesh is
depicted in Figures 13, 14 and IS: the improvements due to residual averaging and the multigrid
scheme are shown in the latter two figures.
Figures 8 and 9 provide a comparison of overall
convergence rates of the triangle code and FL052
on 128x32 meshes. The average residual reduction
per mltigrid cycle of 0.9272 for the triangle
code is somewhat slower than that achieved on
quadrilaterals. However, it m a t be recalled that
the triangular element mesh has twice as many
cells as the quadrilateral mesh.
a
Time stepping on all grids uses the same
Courant number along with a variable time step and
identical residual averaging. Enthalpy damping
is confined to the finest grid, while the dissipative operator D is mdified fos economy on the
coarse grids to a first order accurate form, which
simply uses 2nd differences with a fixed coefficient.
8.
Thus the calculations verify that on equivalent meshes solutions can be achieved with equivalent accuracy and convergence rates using either
triangular or quadrilateral elements. On a
triangular mesh, discretization errors which
depend on the distance between cell centers arise
when flow variables between t w o centers are taken
8 8 the average of the two values BL rhese centers.
However, discretization errors which depend on the
size of the cell face also arise due to approximations when integrating around the boundary of the
cell. If a given area is covered by N quadrilaterals with cell spacing h, then 4 / 3 6 N equilateral triangles are required to cover the same
area, if an equal spacing h is maintained between
the cell centers. Om the other hand, if we retain
the same side length h, then 4 1 6 N equilateral
triangles are required. Thus we cannot precisely
characterize a triangular mesh of equivalent
accuracy to a given quadrilateral mesh. The
numerical calculations used twice as many elements
in the triangular mesh. In subcritical cases the
triangle scheme also realized greater accuracy, as
measured by the drag coefficient, suggesting that
the size of the triangles could have been
increased for equivalent accuracy. Also the
amount of computation is reduced by the need to
calculate fluxes on 3 rather than 4 sides per
cell.
Results
A typical triangular element 0-mesh about an
airfoil is shown in Figure 4. The grid contains
128 nodes around the airfoil and 32 in the radial
direction. We refer LO it as B 128x32 grid,
although there are actually 256x32 triangular elements. The equivalent quadrilateral 0-mesh with
the same number of mesh nodes is depicted in
Figure 5. Figure 6 shows the calculated pressure
distribution about a NACA 0012 airfoil at Oo
angle of attack and Mach number 0.8.
The artificial dissipation is successful in eliminating any
preshock oscillations. and the calculated value of
the drag Coefficient of 0.0085 agrees almost
exactly with that obtained from a calculation
using FL052, an Euler solver on an equivalent
quadrilateral mesh, (see Figure 7). Note that
this value also agrees closely with the eaperimentally observed increase in drag (1.e. net wave
drag) at this Mach number over the subcritical
value for CD (mainly viscous drag).
A supercritical lifting case is shown in
Figure 10 where the NAWL 0012 airfoil has been
given an angle of attack of 1.25O at the same Mach
number.
Figure 11 shows the calculated pressure
distribution about a KORN airfoil at Oo angle of
~
5
9.
Conclusions and Further Work
We conclude that for a given desired accuracy
numerical solution to the Euler equations in 2
dimensions can be obtained using a regular
triangular mesh vith roughly the same amount of
computing effort as vith a quadrilateral mesh.
This indicates that triangular grids m y provide a
useful tool for calculating flows past complex
configurations.
a
b,
The existing code is presently being rewritten
in a more general form, using an indirect
addressing system to locate cell centers and their
corresponding mesh nodes, in order to allow the
use of completely irregular meshes. Another version where the flow variables are stored at the
cell corners is also under development. Planned
extensions of this work include patching of
triangular grids with quadrilateral grids, and
extension to three dimensional grids with triangular prisms or tetrahedra.
Figure 1
Numbering Scheme for a Regular Triangular Grid
References
1.
Jameson. A., "Solution of the Euler Equarions
by a nultigrid Method", Applied Mathematics
and Computation, 13, 1983 pp. 327-356.
2.
Jameson, A., Schmidt, W., and Turkel, E.,
"Numerical Solution of the Euler Equations by
Finite Volume Methods Using Runge-Kutta Time
Stepping Schemes", A I M Paper 81-1259, 1981.
3.
Jameson, AI, and Baker, T.J., "Solution of
the Euler Equations for Complex
Confinurations". Proc. A I M 6th Comoutational
Fluid Dynamics Conference, Danvers, 1983, pp.
-
.~
293-302.
4.
Berger, K . J . and Jameson, A . , "Automatic
Adaptive Grid Refinement for the Euler
Equations", Princeton University, MAE Report
No. 1633, October 1983.
5.
Vermeland, R.E., "Solution of the
Two-Dimensional Euler Equations on a Hybrid
Mesh", Princeton UniVeKSiCy, Department of
Mechanical and Aerospace Engineering Report
No. 1679, August 1984.
6.
P e l r , R.B.,
7.
"Transonic Flov Calculations
Using Finite Elements", Ph.D. Thesis,
Princeton University 1983.
COORDINATE
URECTION 3
IbLIGNED)
Moravetz, C.S, "On the Non-Existence of
Continuous Transonic Flows Past Profiles",
Corn. Pure Applied Math, 9, 1956, pp. 45-48.
Figure 2
Indentification of the 3 Coordinare Directions
(2)
for Calculating the Pressure Switch Eik
6
Figure 3
I n t e r p o l s t i o n of Coarse G r i d R e s i d u a l s t o
the Fine Grid
Figure 4
Figure 5
128x32 Regular Triangular Grid about a
NACA 0012 A i r f o i l
Equivalent 128x32 Quadrilateral Grid
I
0
?
e-
...,.-.._
..
.
.
I.
,..
- .
d-
N4Cii
"#in
0012
0.800
RLP"R
CL
0.0
co
GRIO
128x32
NCIC
H1CR OCl2
nacn
0.500
CL
-0.0000
G R i C iZ5X32
0.0
o.ooes
in
0.0
200
RLi"Fi
cc
0.0
0.0085
NCIC
C"
C.oL0:
200
Figure 6
Figure 7
C a l c u l a t e d Pressure D i s t r i b u t i o n
U s i n g T r i a n g l e Code
C a l c u l a t e d Pressure D i s t r i b u r i o n Using FLOIZ
NRCR 0012
*9:"
0.800
R E S I D > 0.2530+01
"OR*
193.00
WPHR
NRCII 0012
Mac"
0.0
RETE
0.800
RES101 0.1790101
"OR"
193.00
RES102 0.7VSD-06
0.9212
Figure 8
0.0
RES102 0.1250-07
RILPIIR
RRTE
0.3099
Figure 9
Convergence Rate of n o 5 2
Convergence Rate of K u l t i g r i d T r i a n g l e Code as
U e a s u r e d by t h e RHS Value of a d a t ( d e n s i t y r e s i d u a l )
and Build-up of S u p e r s o n i c P o i n t s .
8
*.......
_..-.+
_
_
.
.
e
+
/-.
.
.. -.
-..
a.
1
.
-a.
.. 2.
-..2..
.....
..:.
1.
z.
--i
d
D
0.800
LL
GRID
0.3522
I28332
RLPHR
1.25C
LO
0.022:
NCTC
201
CK
t
KORN R I R F O I L
NRCR 0012
MRC*
I'
Hac*
CL
-0.0373
GRID
0.750
0.62u5
128x32
0.0
RLPHR
CO
NCIC
0.0003
Figure 10
Figure 1 1
Supercritical Lifting case
KOEN Airfoil
32x8
cv
-0.!~5$
800
a4170
OOOM
-a0047
0.4244
aooig
-0.0040
64x16
0.4272
0.0009
-0.0039
96x24
04290
O.WO3
-0.0039
128x32
0.4296
0.MX)I
-0.0039
Table 1
.. .
HOLM
CL
GRID
__ ..
0.500
0.U296
128132
Calculated Force Coefficients as a Function of Grid
Size for the Subcritical Lifting Case
RLPHR
CD
3.000
0.0001
NCTC
Cr
-S.ocjs
300
Figure 12
Subcritical Lifting Case
9
D
7
NFCR 0012
natH
0.8OC
R E S I D 1 0.2630.01
YORn
798.00
GRID
LVX16
NU:e
RLPW
ncCH
0.t
R E S 1 0 2 O.20:C-C3
RRIE
C.58SI
Figure 13
Y"i%
796.00
GRiS
SYXI6
nmfi
0.0
FESIG2 0.278D-OE
RPlE
O.56CE
Convergence Rate with Residual Smoothing
NRCR 0012
6.0
R E S 1 0 2 0.3U-3-CY
DLPHC
RRiE
0.800
Figure 1 4
Convergence Rate as Measured by the RMS Value
of a p / a t ( d e n s i t y r e s i d u a l ) and Build-up of
Supersonic P o i n t s f o r the Standard Triangle Code
IlRC"
0.800
R E S I D 1 O.iY3D*Ol
YORK
199.00
GRID
6YX16
OOl2
R E S I D ) o.lYSD10l
0.85.:
Figure 1 5
Convergence Rate w i t h Multigrid
10