A Cell-Vertex Multigrid Method For The Navier-Stok

Download as pdf or txt
Download as pdf or txt
You are on page 1of 42

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/24296656

A cell-vertex multigrid method for the Navier-Stokes equations

Article · February 1989


Source: NTRS

CITATIONS READS

41 157

1 author:

Rolf Radespiel
Technische Universität Braunschweig
391 PUBLICATIONS 3,823 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

SFB 880 - Fundamentals of High Lift for Future Civil Aircraft View project

Metropolitan Aircraft - Bürgernahes Flugzeug View project

All content following this page was uploaded by Rolf Radespiel on 09 June 2014.

The user has requested enhancement of the downloaded file.


NASA Technical Memorandum 101557

A CELL-VERTEX M U L T I G R I D METHOD FOR THE


N A V I ER-STOKES EQUATIONS

(hASA-TH-lOlS57) A CkLL-YEb2EX LIOL'LIGBID N89-201Cl


L E I P G C FOR 9 B E E A P I E R - 5 I G K E S I C t i A ' I I G H S
OASA) 41 CSCL O l A
Unclas
G3/02 CISESt3

R. Radespiel

JANUARY 1989

National Aeronautics and


Space Administration
Langley Research Center
Hampton,Virginia 23665
I

Summary
A cell-vertex scheme for the Navier-Stokes equations, which is based on central difference
approximations and Runge-Kutta time stepping, is described. Using local time stepping,
implicit residual smoothing, a multigrid method, and carefully controlled artificial dissipative
terms, very good convergence rates are obtained for a wide range of two- and three-
dimensional flows over airfoils and wings. The accuracy of the code is examined by grid
refinement studies and comparison with experimental data. For an accurate prediction of tur-
bulent flows with strong separations, a modified version of the nonequilibrium turbulence
model of Johnson and King is introduced, which is well suited for an implementation into
three-dimensional Navier-Stokes codes. It is shown that the solutions for three-dimensional
flows with strong separations can be dramatically improved, when a nonequilibrium model of
turbulence is used.

Introduction
In recent years, considerable advancement has been achieved in the numerical solution of
the Euler and Navier-Stokes Equations. Nevertheless, existing computer codes for the solution
of the three-dimensional Navier-Stokes equations require large computing times. In order to
make Navier-Stokes solutions useful in the flight vehicle design process, substantial improve-
ments in efficiency and accuracy of the algorithms have to be made.
Finite-volume methods based on explicit Runge-Kutta time stepping schemes have been
shown to be very efficient in the solution of the Euler equations governing inviscid flows
[1,2]. They are in widespread use now. More recently, they have been extended for a solution
of the mass-averaged Navier-Stokes equations. The convergence of the basic scheme to the
desired steady state solution, however, slows down considerably because of the time step lim-
itation associated with the small mesh cells which are necessary to resolve the thin shear
layers. This drawback of the explicit scheme can be overcome by applying several accelera-
tion techniques, namely local time stepping, implicit residual averaging and a multigrid algo-
rithm. Successful applications of these techniques have been reported in [3,4] for two-
dimensional flows and there has been a partial success in three dimensions also [5].
In the present work a recently developed efficient and robust finite-volume multigrid
scheme for the two- and three-dimensional Navier-Stokes equation is presented. The scheme
employs a cell-vertex discretization rather than the usual cell-centered discretization because
. an analysis of the discretization errors indicates lower errors for meshes with high stretching
or discontinuities in slope. The accuracy of the treatment of the far-field boundary is
improved by representing the lifting surface by a vortex and using the induced flow field to
determine the free-stream conditions. As the scheme is based on central differencing, an
artificial viscosity model is used, which is particularly designed for high aspect-ratio cells. The
boundary treatment of the dissipative terms, which is important for an accurate prediction of
the skin friction, is discussed. A multistage Runge-Kutta time-stepping algorithm is used to
advance the solution in time. The acceleration techniques, which are applied to obtain faster
-2-

convergence, are discussed. For an accurate prediction of turbulent flows with strong separa-
tions, a modified version of the nonequilibrium turbulence model of Johnson and King [6] is
introduced, which is well suited for an implementation into three-dimensional Navier-Stokes
codes.
Numerical results for subsonic, transonic and supersonic flows around airfoils and finite
span wings are given. The accuracy and the convergence behavior of the scheme are investi-
gated by varying the grid density and the coefficient of the artificial dissipation. The influence
of the turbulence model on the numerical results is shown for several flow cases. It is
demonstrated that the solutions converge very rapidly to the steady state. Typically, engineer-
ing accuracy is obtained within 40 multigrid cycles for high Reynolds number transonic airfoil
flows on a 321x65 mesh. For finite-span wings, converged solutions are obtained within 50-80
iterations on a mesh with 288x65~49points.

half span
speed of sound
skin-friction coefficient
pressure coefficient
drag coefficient
lift coefficient
Courant number of time stepping scheme
Distance from the wall
total internal energy
unit vectors of Cartesian coordinate system
coefficient of the second difference dissipation
coefficient of the fourth difference dissipation
Mach number
unit vector of outer normal
pressure
Prandtl number
cell face vectors in the 6-, q-,and [-directions
time
temperature
Cartesian velocity components
.
volume
vector of conserved quantities = (p pu pv pw P E ) ~
Cartesian coordinates
angle of attack
ratio of specific heats, intermittency factor
-3-

central difference operator, boundary layer thickness


coefficients of implicit smoothing according to (26)
v. Karman's constant = 0.4
spectral radii of the Jacobian matrices associated with the {-,q-,[-directions
viscosity according to (4)
density
nonequilibrium factor in J.-K.-mod. turbulence model
shear stress
exponent in (12) and (27)
magnitude of vorticity
curvilinear coordinates

discrete quantity
location of maximum of function g, equation (36),across boundary layer
location of maximum of function F, equation (34), across boundary layer
turbulent
value at the wall
value infinitely far away from the aerodynamic surface

Governing Equations
The integral form of the mass-averaged Navier-Stokes equations using nondimensional
variables can be written as

? +JJF3dS = 0
~ J J J VdV
at v av
where
3 = (p pu pv pw P E ) ~ .
is the vector of conserved quantities with p, u, v, w, E denoting the density, the Cartesian
velocity components and the specific total internal energy. The quantity V denotes an arbi-
trary control volume with the boundary aV and the outer normal 3.The flux tensor F may be
divided into its convective part F, and its viscous part F, as
- - -
F = F, - F,
I \

pu +pv; +pwx
(pu"+p)Z + puv ly + puw +1,
7)

-
Fc = puvx + (pv2+pG + pvw 1, 7+

puw; + pvw +ly + (pw2+p)j


~ (puE+upC + (pvE+vp)< + (pwE+wp)X a
-4-

-
F, =

The equation of state for an ideal gas is used to calculate the pressure, p, and the temperature,
T. That is,
2 2
P = (Y-l)p( E - (U +V 112 1 , T = PIP (2)
where y denotes the ratio of specific heats. The elements of the shear-stress tensor and the
heat-flux vector are given by the constitutive equations for a Newtonian fluid as follows:
ox, = 2pux - 2/3p(ux+vy+wz,
on = 2pvy - 2/3p(ux+vy+xz)
a,, = 2pv, - 2/3p(u,+vy+x,)
oxy = a y x = WY+VX) (3)
oxz
-
- a z x = p(u,+w,>
by2 = a z y = p(vz+wy)
9, = -kaT/ax, qy = -kaT/ay, qz = -kaT/& .
The nondimensional viscosity, p, is assumed to follow an empirical power law

and the heat conductivity is

For turbulent flows, the laminar viscosity, p, is replaced by p+pt and p/Pr is replaced by
p/Pr+pt/Pr,, where the eddy viscosity, pt, and the turbulent Prandtl number, Prt, are provided
by a turbulence model. The turbulence models used in the present work are described in a
subsequent paragraph.

Spatial Discretization
The discretization of (1) follows the method of lines, i.e. discretization in space and time
are done separately. The domain around the aerodynamic body is partitioned with hex-
ahedrons. The discrete values of the flow quantities are located at the vertices of the mesh
-5-

cells. The eight cells surrounding a vertex form a super cell as shown in Fig. 1. The integral
equation (1) is approximated by the spatial discretization yielding

In the following, the net inflow of mass, momentum and energy associated with the convec-
tive terms Gci,jk, the viscous terms dk and the artificial dissipative terms B i j , k are dis-
cussed. For this purpose the super cell around the point (ij,k) for the simplified case of a
plane mesh is shown in Fig. 2.
The surface integral of (1) is evaluated for each component cell using an arithmetic aver-
age of the flux quantities at the vertices to determine the values on each of the boundary seg-
ments. Then the resultant convective inflow of mass, momentum and energy associated with
point (ij,k) is computed by summing the contributions of the component cells. According to
[7], the scheme is at least first-order accurate, if two conditions are satisfied. (1) The distribu-
tion of the normal vector over each surface segment has to be smooth, that is, if the grid is
refined, the four grid points defining a segment tend to lie within a plane. (2) The shape of
each surface segment has to approach a parallelogram with grid refinement. Note, that these
conditions still allow discontinuities in the slope of the grid lines and discontinuities in the
rate of grid stretching. On completely smooth meshes, the discretization is second-order accu-
rate.
Next, the viscous fluxes required to determine the solution at the point (i,j,k) are approxi-
mated using the auxiliary cell with the dashed boundary shown in Fig. 2. The viscous fluxes
contain first derivatives of the flow variables, which are computed using a local transformation
from Cartesian coordinates to the curvilinear coordinates c,q,[, i. e.

The derivatives $6, $q and $c are approximated using finite differences, whereas the cell face
vectors $e, 4,sQ which are written as

and the volume V are used to compute metric derivatives. If QX is to be approximated at


(i+l/2,j,k) we obtain

where 9. Sq, 66 denote central difference operators in the curvilinear coordinate directions.
In practice, coordinate grids for the resolution of viscous shear layers are highly stretched in
-6-

the direction normal to the layer and therefore, an one-dimensional error analysis of the
discrete approximation to the viscous terns can give useful insight. In [4], it was shown that
the present scheme is first-order accurate on general stretched meshes and second-order accu-
racy is obtained on smoothly stretched meshes. In the present version of the code, the viscous
terns have been simplified by taking into account gradients in the direction normal to the
viscous shear layers only (thin layer approximation).
In order to prevent odd-even point decoupling and oscillations near shock waves,
artificial dissipative terms are added to the governing discrete equations. The artificial dissi-
pation model considered in this paper is based on the work of Jameson, Schmidt and Turkel
[l]. A blend of fourth and second differences is used to provide third-order background dissi-
pation in smooth regions of the flow and first-order dissipation at shock waves, and is given
by

i The second and fourth difference operators read

where A<, Vg are forward and backward difference operators in the 6-direction. In order to
I
avoid excessively large dissipation levels for cells with high aspect ratios and to maintain
good damping properties of the scheme, a variable scaling factor of the dissipative terms is
employed, which is a function of the spectral radii of the Jacobian matrices associated with
the c,q,c directions and accounts for varying cell aspect ratio
-
, kt i+lRJ.k = kg all2,j.k * @i+1/2,j.k 9 (12)

where q + l Q j , k denotes the velocity vector and c is the speed of sound. The use of the max-
imum function in the definition of Q, is important for grids where h,l/he and k& are very
large and of same order of magnitude. In this case, if these ratios are summed rather than tak-,
ing the maximum, too large dissipative terms are obtained, which may destroy the conver-
gence of the time-stepping scheme. In general, coordinate grids around three-dimensional
configurations exhibit larger variations of the cell aspect ratio than two-dimensional grids.
Therefore, we have observed somewhat less freedom in the choice of the exponent o in three
dimensions. It has been found that the choice ~ 0 . yields5 a robust three-dimensional scheme.
-7-

The coefficients E(') and d4)use the pressure as a sensor for shocks in the flow field. They are
defined as

Pi-lj,k - 2pij,k Pi+lj,k


I
"ij,k = I
Pi-1j.k + 2pi,j,k Pi+lj.k

where k(2) and k(4) are constants. The dissipation operators in the q- and c-directions are
defined in a similar manner.

Boundary Conditions
At subsonic inflow/outflow boundaries, locally one-dimensional flow normal to the boun-
dary is assumed. According to [8], the governing flow equations are linearized around the
values at the end of the previous time step, and the characteristic variables corresponding to
outgoing characteristics are extrapolated from the interior. The characteristic variables
corresponding to incoming characteristics are determined from the free stream.
For two-dimensional flows around airfoils, it is well known that the free-stream condi-
tions have to be specified by taking into account the circulation around the airfoil to obtain
correct values of lift and drag for a given angle of attack. The influence of an improper
determination of the free stream at the far field may be felt even if the distance to the far field
is as large as 20-50 chord lengths. The usual practice in calculating the free-stream condi-
tions for flows around finite-span wings has been to use simply the undisturbed onset flow
conditions. In fact, if the far field is sufficiently far away from the wing surface, the distur-
bances created by the wing behave like those coming from a point singularity rather than a
line singularity as in two-dimensional flow. Hence, a disturbance would decay more rapidly
in three-dimensional flow than in two-dimensional flow, thus allowing the use of a smaller
distance to the far-field.
More accurate results can be expected if the free-stream conditions used in three-
dimensional flow computations take into account disturbances created by the aerodynamic sur-
faces. Klunker [9] shows that the small-disturbance potential at the far field is dominated by
terms representing the lift of the aerodynamic surface. He derives an equation in which the
potential at the far field is determined by integrating the contributions of singularities
representing the local circulation over the entire wing surface. This approach has been used in
[ 101 to determine the free-stream conditions in an Euler code.
If free-stream conditions are to be specified for general configurations a numerical
integration over the aerodynamic surface for each point in the far field seems to be untract-
able. In particular, if multiblock structured meshes are to be used, this integration requires
considerable amounts of programming effort and storage. The determination of the free
stream is much simpler if the wing is replaced by a horse-shoe vortex. The components of the
-8-

disturbance velocity vector induced by a horse-shoe vortex in compressible flow are*

U =

v=-[
-r z+a
(1+-)
X
- z-a (1+-)
x
+ p2x z+a z-a
2~ (z+a)2+y2 .Iw, (z-a>2+y2 .Iw- x2+p2y2K

where
P=.\11-Mm2,
y+= x2+p2(z+al2 + p2y2 ,
y- = x2+p2(z-al2 + p2y2 .
The circulation of the vortex is related to the lift of the wing by

the quantity a denotes the half span of the horse-shoe vortex and ,c, 1s the mean chord of
the wing. The position of the horse-shoe vortex with respect to the Cartesian coordinates x,y,z
is shown in Fig. 3. The induced velocity is infinite downstream of the wing where the vortex
crosses the boundary. In reality, however, this is by no means a special point in the flow field.
Furthermore, numerical experience indicates that the treatment of the outflow boundary down-
stream of the wing has no significant influence on the results. Therefore, the following numer-
ical treatment has been employed. The denominators in (15) and (16) which are responsible
for the singular behavior of v and w at downstream boundaries, are limited in their value for
x>O, i.e. (z+a)2+y2 and (z-a)2+y2 are kept at values larger than a/2 for x>O. This means that
the induced velocities are artificially reduced in their value within the distance a/2 around the
vortex line which crosses the downstream far-field boundary.
For all grid points at the far field, the geometric relations in the formulas given above
can be evaluated once and stored. At each time step, these influence coefficients are multi-
plied by the total lift of the wing and added to the onset flow. Thus, the work to obtain the
free-stream velocities is negligible. The pressure and the density of the free stream are then
obtained using the assumption of isentropic flow and constant total enthalpy. The numerical
results obtained with the present determination of the free-stream values show a reduced sen-
sitivity to the far-field distance (see section on numerical results).

* These relations have been dcrivcd in unpublished ’Notes on Linearized Subsonic Wing
Thcory’ by E.B. Klunkcr and K.C. Harder.
-9-

At solid walls, the no-slip condition is enforced by setting

The continuity and energy equations are solved at the grid points lying on the surface assum-
ing an adiabatic wall. For this purpose, a control volume is formed with the four nearest cells
outside and their mirror images inside. When updating the convective terms, the flow vari-
ables at the computational points inside are obtained as the symmetric images of the values
outside. For the viscous terms, the density and the total energy inside are obtained as sym-
metric images, whereas the velocity components are taken as the antisymmetric images of the
values just outside.
It has been found that the computed velocity distributions near the wall may be
significantly affected by the artificial dissipative terms if the stencil of the dissipative terms
normal to the wall is not properly defined. The reasons for this are the high grid stretching
normal to the wall and the steep velocity gradients of the turbulent boundary layer. When
solving the momentum equations at the grid point just outside the wall at j=3, the dissipation
operator needs the dependent variables just inside. The extrapolation formula

rather than the usual linear form

results in much less sensitivity of the solution to artificial viscosity and has had no significant
drawback on convergence in the computations made so far.

Time Stepping Scheme


The system of ordinary differential equations which is obtained by the discretization in
space is advanced in time with a five-stage Runge-Kutta scheme. A hybrid scheme is used
where the physical viscous terms are computed on the first stage and frozen for the remaining
stages and the artificial dissipative terms are evaluated on the first, third and fifth stages of the
scheme. At the (m+l)st stage,

1 n=O J

where qCo) is the solution at the old time level, and At is the time step. The subscripts c and
v refer to convection and physical diffusion contributions to the discrete flow equations. The
stage coefficients, a, are taken as
- 10-

The weighting factors of the artificial dissipation satisfy the condition


CYm=1
They are

Y43 =0 ? Y44 = 7.5 9

Y3
where = 0.56,T. = 0.44. It has been shown in [3], that this scheme has a particularly large
parabolic stability limit. The steady-state solution is independent of the time step and there-
fore, the scheme is amenable to convergence acceleration techniques.
Three methods are employed to accelerate convergence of the basic explicit time step-
ping scheme. These techniques are as follows: (1) local time stepping, (2) implicit residual
smoothing and (3) multigrid.
With local time stepping, the solution at each mesh point is advanced at the maximum At
allowed by stability. Both convection and diffusion stability limits are included in the compu-
tation of At as

where
4 1 P+Pt
Vdiff = max(-, ) -,
3 p r P
and the coefficient Cdiff = 2 is used.
The stability range of the basic time stepping scheme can be extended using implicit
smoothing of the residuals [ 113. For three-dimensional flows the residual smoothing is applied
in the form

where R, is the explicit residual at the Runge-Kutta stage m in (21), and is the final resi-
dual at stage m after the sequence of smoothing in the 5, q and 6 directions.
The use of constant coefficients in the implicit treatment has proven to be satisfactory
(extending the CFL number by a factor of two to three) even for meshes with cells of high
aspect ratio, provided additional support such as enthalpy damping [l] is introduced. How-
ever, the use of enthalpy damping, which assumes constant total enthalpy throughout the flow
- 11 -

field, precludes the solution of problems with heat transfer effects. The need for enthalpy
damping can be eliminated by using variable coefficients Q,, E, and E< that account for the
variation of the cell aspect ratio.
For this purpose, consider a cell where the edge lengths in the 5- and c- directions are
much longer than that in the 7-direction. The explicit time step is limited by the characteris-
tic wave speed in the direction of the short cell edge. It is obvious that, for an extension of
stability, implicit smoothing is required in the 7-direction. If the same implicit residual
smoothing is also applied in the other directions, where the characteristic wave speeds are
much smaller than the stability limit, the damping behavior of the scheme, which is optimal at
wave speeds near the stability limit, is impaired. To overcome this problem, Martinelli [3]
has given a formula, in which the smoothing coefficient is a function of the characteristic
wave speeds. Martinelli's form and variable smoothing coefficients in general are discussed in
detail in the unpublished work of Swanson. Here, it is extended to three dimensions, i.e. q
may be defined as

Finally, a multigrid method is employed which is based on the work of Jameson [2]. For
the multigrid process, coarser meshes are obtained by eliminating every other mesh line in
each coordinate direction. The solution is transferred to coarser meshes by injection. Residu-
als are transferred from fine to coarse meshes by a weighted average over the fine mesh grid
points which are nearest to the point on the coarse mesh. A forcing function is constructed so
that the solution on a coarse mesh is driven by the residuals collected on the next finer mesh.
This procedure is repeated on each succeeding coarser mesh until the coarsest mesh is
reached. Then, the corrections are transferred to the next finer mesh by trilinear interpolation.
A fixed W-type cycle is used to execute the multigrid strategy. The robustness of the mul-
tigrid scheme is greatly improved by smoothing the total corrections before the solution is
updated. That is,

where

The quantity A$F is the solution correction from the finest grid, and A$c is the resultant
solution correction from the coarse grids. The smoothing reduces high frequency oscillations
introduced by the linear interpolation of the coarse mesh corrections and hence, convergence
of the scheme is obtained for a much broader range of dissipation coefficients. The factored
scheme of (26) with constant coefficients (E{ = e,., = EC = 0.2) is used for this smoothing.
Additional robustness and efficiency of the multigrid scheme is achieved by computing more
than a single time step on the coarse meshes. In the first place, the use of multiple time steps
on the coarse meshes improves the convergence, because the solution is further advanced in
- 12-

time on the coarse meshes. Secondly, multiple time stepping improves the damping charac-
teristics on the coarse meshes. Thus, smaller values of coarse-mesh dissipation may be used
which results in an improvement of the convergence to the steady state. Finally, the applica-
tion of a Full Multigrid method provides a well conditioned starting solution for the finest
mesh being considered.

Turbulence Models
In the present work two turbulence models are considered. The first model is the widely
used algebraic model of Baldwin and Lomax [12], which is based on the eddy viscosity
hypothesis. The eddy viscosity is defined as
Pt = fin(Pt,i, ~t,o) - (29)
In the inner region of the boundary layer,
Pzi = P(K D d)* R (30)
where K is von Karman's constant, f2 is the magnitude of the vorticity, and the coordinate d
is the distance from the wall. The quantity D is the van Driest damping factor

V tJw Ll,max
D = 1 - exp(-y+/26l -,>
v+ = d
,
PW

Note, that the maximum laminar shear stress across the layer, z ~ rather
, than
~ the~ shear
~
stress at the wall is used to define the damping term in order to prevent vanishing eddy
viscosity at the separation points. In the outer region,
Pt,, = 1.6 * 0.0168~F,*~y,
where

Fmaxis the maximum value of


F=DRd (34)
across the layer, and is the value of d at which it occurs. The quantity Udif is the
difference of the flow speed across the shear layer. The intermittency factor y is
Y =
1
(35)
1+5.5 (0.3-d/d,,)6
This turbulence model will be denoted by B.-L..
Frequently, it has been observed that in separated flows, the assumption of an equili-
brium of the turbulence is not valid. Both convection and diffusion of turbulence have to be
taken into account to allow for a more accurate determination of the turbulent stresses. For
this purpose, Johnson and King [6] have devised an extended eddy viscosity model for
separated wall boundary layers. They postulate that the viscosity of the inner layer near the
- 13 -

wall and the shape of the eddy-viscosity distribution in the outer layer may be accurately
described with algebraic relations assuming equilibrium. The equilibrium distribution of the
outer layer, however, is multiplied with a nonequilibrium factor in such a way that a transport
equation for the maximum shear stress in the layer is fulfilled. Thus, the complete model is a
nonequilibrium model of turbulence. In the following, a slightly modified version of this
model is described. The model is particularly convenient for use in three-dimensional Navier-
Stokes codes by using the Baldwin-Lomax formulation for the outer viscosity as proposed in
r133.
First, instead of using the maximum turbulent shear stress as originally proposed in [6],a
turbulence reference quantity,
g, = (Q cLt~P)m (36)
is introduced, where the index m denotes the maximum of g across the shear layer. Note, that
g, is invariant to the coordinate system used, and therefore, it can be easily used in Navier-
Stokes codes. According to 163, the eddy viscosity is defined by an exponential blending of
the inner and the outer viscosity
.

The inner viscosity is written as

L J

The use of the maximum function in (38) is important for numerical stability on coarse
meshes, where the turbulent boundary layer is not well resolved. According to [14], a value of
17 has been chosen for A'. In the outer layer, the original model of Johnson and King uses
the flow velocity at the edge of the boundary layer and the displacement thickness to define
the outer viscosity. The displacement thickness, however, is not easily determined in three-
dimensional flows. In order to avoid numerical problems associated with searching through
the boundary layer in three-dimensional flows, the viscosity in the outer layer is determined
by the relation of Baldwin and Lomax
vs0= CJ 1.6 * 0.0168pFW*,y. (39)
The quantity CJ is the aforementioned nonequilibrium factor, and Fw&, is defined according to
(33). The intermittency factor, however, is redefined as
.
I
1
Y=
1+5.5(d/6)6
- 14-

where 6 denotes the thickness of the boundary layer. In [15] it has been demonstrated that
~

6 = 1.9 d,
is a good approximation for wall boundary layers, and this relation is incorporated into the
model. Furthermore, the behavior of the model in separated flow can be improved, if the dis-
tance d in (34) for the function F is computed according to Fig. 4 rather than taking the dis-
tance from the wall.
Assuming that g is proportional to the kinetic energy of turbulence, the distribution of g
,
over the aerodynamic surface is determined by a differential equation given in [6]. For three-
dimensional flows, the equation may be written as

Here, u,, v,, and w, are the Cartesian velocity components at the location of g,. The index
eq,m denotes the equilibrium value at the location of ,,g and the length scale, L,,,,
is defined
I
as
I
L, = 0.44, for d,,,/6 I 0.225 , (43)
L, = 0.096 for d,,,/6 > 0.225 .
~

The diffusion term, D, is defined as


a2 max(0, 0 l ’ ~ - 1 )
D= (44)
6(0.7 - (d/6),) *

The model constants al=0.25 and a2=0.5 are used here. Equation (44) states that the diffusion
of turbulence is neglected in retarded flows, where 04,whereas in recovering regions, the
diffusion is assumed to be proportional to d n - l . If g
, is replaced with grn=gm-ln,a linear
equation is obtained as

(45) is valid along a surface which is determined by the location of the maximum of the refer-
ence quantity g. For numerical convenience, (45) is solved on the surface grid of the aero-
I dynamic body rather at the location of .g, The misalignment between these surfaces creates
errors in the convective terms of (45). These errors are reduced by replacing the velocity com-
ponents u,, v,, and w, in (45) by their projection onto the body surface. It is believed, that
the remaining errors, which are due to the skewness of the grid lines crossing the shear layer,
are small in comparison with the uncertainties present in the basic assumptions of this tur-
bulence model. Note, that the time derivative in (45) is retained, so that a time stepping
scheme can be used for the numerical solution. Using Gauss’s theorem, a cell centered
finite-volume discretization with central differencing is applied to (45). Because g, does not
vary normal to the surface, only a layer of one cell thickness around the surface is required to
- 15 -

solve (45). Fourth difference dissipative terms are used to damp the discrete equations. The
explicit five-stage Runge-Kutta scheme (21) is used for time stepping. The linear source term
in (45) is treated implicitly, and thus, CFL-numbers around 3 can be used. The convergence
to the steady state is further enhanced by local time stepping. For the computation of flow
cases with strong separation the turbulent viscosity should be rapidly adjusting to the flow
field in order to allow for convergence. Therefore, 10 explicit time steps are executed at each
time, when the turbulence model is updated. The computational time for the solution of (45)
is very small in comparison with the total effort, because only a single equation for a quantity,
which does not change normal to the surface, is involved.
Once values of g, have been computed by the numerical solution of (43, they can be
compared to values of g, from (36). The ratio of both these values is then used to update 0,
as described in [131; and eventually, convergence is obtained. The nonequilibrium turbulence
model described herein will be denoted with J.-K.-mod. .

Numerical Results
The cell-vertex multigrid scheme described here has been applied to a wide variety of
two- and three-dimensional flows. The two-dimensional subsonic laminar flow around NACA
0012 airfoil at the free-stream Mach number, M, = 0.5, the angle of attack, a = Oo, and the
Reynolds number based on the chord length Re, = 5000 is considered first. A C-type coordi-
nate mesh of 257x65 points is used, where 129 points are distributed over the airfoil surface.
The first spacing away from the wall is 5x104 chord lengths. Fig. 5 shows the convergence
history, the streamlines, and the distributions of pressure and skin friction. The computed
values of pressure drag and skin-friction drag and the separation point agree well with the
results of Swanson and Turkel [16]. The solution converges rapidly, and, noting additional
convergence criteria given in Table 1, a sufficiently accurate solution is obtained in less than
40 multigrid cycles on the fine mesh.
The transonic turbulent flow over the RAE 2822 airfoil is used to demonstrate the con-
vergence behavior and the accuracy of the method for high Reynolds number turbulent flows.
A mesh of 385x65 points is used, where 257 points are distributed over the airfoil surface,
and the first spacing away from the wall is chord lengths. Frequently, a mesh with
321x65 points has been used also, which has the same distribution of grid points at the sur-
face, but less points along the wake line. In Fig. 6, the convergence behavior is displayed for
M, = 0.73, a = 2.79', Re, = 6 . 5 ~ 1 0 ~Fig.
. 6a shows the influence of the fourth-difference
dissipation on the convergence of the solution. The robustness of the scheme is demonstrated
by the fact that the dissipation coefficient is varied by an order of magnitude without destroy-
ing the convergence of the scheme. Note, that the results of Fig. 6a have been obtained using
single time steps on the coarse meshes. The effect of multiple time steps on the coarse
meshes is shown in Fig. 6b. If two time steps are used, the convergence is significantly
improved, and a residual reduction of 10 orders of magnitude is obtained within 340 multigrid
cycles on the fine mesh. The convergence criteria in Table 1 indicate that sufficiently
- 16-

converged solutions are obtained within 40 multigrid cycles. Fig. fx shows the convergence
with the nonequilibrium model of turbulences, which is slightly slower than with the B.-L.
model.
Next, the accuracy of the scheme is examined using a variation of the grid density. For
this purpose, a coarser mesh of 193x33 points is created by omitting every other point in each
coordinate direction. Additionally, a very fine mesh of 577x97 points has been used. The
variation of the coefficients for lift, pressure drag, and friction drag with number of mesh
points N is presented in Fig. 7. The effect of the second and fourth difference dissipation is
also indicated. On coarse meshes, the discretization error is obviously dominated by artificial
dissipation. The finer meshes allow the extrapolation of the coefficients to their values for an
infinitely fine mesh. For the 385x65 mesh (N-'-4.0~10-~), the predicted lift is within 1.5 per-
cent, the pressure drag is within 3 counts, and the friction drag is within 0.3 count of the
extrapolated values. For the 577x97 mesh (N-'-1.8~10-~), the predicted lift is within 0.5 per-
cent, the pressure drag is within 1 count, and the friction drag is within 0.1 count of the extra-
polated values. Fig. 8 shows pressure and skin-friction distributions for different grid densi-
ties, with a comparison to the experimental data of [17]. The main features of the flow are
essentially captured with the 193x33 mesh. The differences between the solutions on the
385x65 and 577x97 meshes are small.
In Fig. 9, the effect of the turbulence model is presented. For both the equilibrium B.-L.
model and the nonequilibrium J.-K.-mod. model, the position of the shock and its strength
agree fairly well with the experimental data. The measured lift coefficient, however, is in
better agreement with the results of the J.-K.-mod. model than with B.-L. The skin friction
distribution of the J.-K.-mod. computation shows relatively high values downstream of the
shock, which are not observed experimentally. The high skin friction in the recovery region is
caused by the definition of the inner viscosity, (38), which is scaled by the maximum of the
quantity g rather than taking local quantities as in (30). The definition (38) should be refined
in future work.
In the following, the behavior of the turbulence models is presented for test cases with
strong separation, which have been used recently to compare various computer codes for
viscous transonic airfoils [18]. In Fig. 10, the flow around the RAE 2822 airfoil at
M, = 0.75, a = 2.81°, Re, = 6 . 5 ~ 1 0is~considered. With the €3.-L. model, the shock location
is predicted 8% of the chord length downstream of the experimental data. The agreement of
computation and experiment is considerably improved, when the J.-K.-mod. model is used,
however, there are still some differences remaining. The bump in the pressure distribution
downstream of the shock, as obtained with the J.-K.-mod. model, indicates that the height of
the separated region, as well as the recovery process downstream are somewhat overpredicted.
Similar trends are observed for the next test case, the NACA 0012 airfoil at M, = 0.799,
a = 2 . 2 6 O , Re, = 9 . 0 ~ 1 0 ~Fig.
. 11 shows much better agreement in pressure distributions and
lift and drag values, when the J.-K.-mod. model is used instead of the B.-L.
- 17 -

Drag polars have been computed for the NACA 001 airfoil at M, = 7 and
Re, = 9 . 0 ~ 1 0 ~The
. results of these computations are shown in Fig. 12 (lift versus angle of
attack) and Fig. 13 (lift versus drag). Experimental results of Hams [19] with wind-tunnel
corrections included are also displayed. Good agreement of computations and measurements is
. obtained at lower angles of attack. At high angles, the B.-L. model yields considerably higher
lift values, than the J.-K.-mod. In particular, the J.-K.-mod. model predicts maximum lift with
CL,max = 0.728 at a = 6', whereas the maximum lift is CL,max= 0.87 at a = 7' with the B.-L.
model. There is no experimental value of the maximum lift for this particular Mach number
given in [19]. The experimental values of the maximum lift at Mach numbers below and
above, CLsnax= 0.82 for M, = 0.65 and CL,max= 0.55 for M, = 0.76 indicate that the result
of the J.-K.-mod. model is close to the probable result of the wind tunnel.
The last two-dimensional test case is the supersonic turbulent flow around NACA 0012
airfoil at M, = 2.0, a = 0', Re, = 10.0~10~. The convergence history and the distributions of
pressure and skin friction are shown in Fig. 14. Generally, good convergence rates are
obtained for Mach numbers up to 2. For higher Mach numbers, no solution could be obtained
with the present scheme. In Fig. 15, the Mach contours and the mesh near the airfoil are
shown. No attempt has been made to adapt the mesh to the bow shock, and therefore, the bow
shock is not very well resolved in this computation.

The three-dimensional version of the code has been applied to the flow over the finite-
span ONERA-M6 wing and to the three-dimensional flow over an airfoil mounted in the wind
tunnel, where the viscous layer along the wind-tunnel side wall has been included in the simu-
lation. In the present report, only the computations of transonic flows over the ONERA-M6
wing are presented. The results for the interaction between a wing and a side-wall boundary
layer are given in a subsequent paper.
The computational domain around the wing is discretized using C-type topology in the
streamwise direction and an 0-type topology in the spanwise direction with 289x65~49points.
The distance to the first grid point away from the wall is local chord lengths. The struc-
ture of the mesh is shown in Fig. 16a. There is considerable skewness in the wing-tip region,
as indicated in Fig. 16b. A somewhat coarser mesh has been also used with 193x49~33
points. In the three-dimensional computations done so far, only a single time step is executed
on the coarse meshes. Fig. 17 shows the convergence behavior of the three-dimensional
scheme for M,=O.84 and Re,=ll.x106. Two flow cases, namely a = 3.06' with attached flow,
and a = 6.06' with a strong separation on the upper surface are considered. For the
attached-flow case the influence of the fourth-difference dissipation on convergence is
displayed in Fig. 17a. The convergence is not as good as for the two-dimensional cases. The
largest residuals occur at the wing tip, where the grid is highly skewed. Therefore, a consid-
erable improvement of the convergence can be expected, if the grid quality were improved in
this region. The effect of multiple time steps on the coarse meshes for the attached-flow case
is shown in Fig. 17b. The decrease of residuals versus multigrid cycles is not much affected
- 18 -

by mu iple time stepping on the coarse meshes. Apparently, talereare hig frequency oscilla-
tions on the fine mesh, which are only slowly damped due to the grid skewness, and this
behavior is not improved using more work on the coarse meshes. The convergence of the glo-
bal flow field, however, is improved as shown in Table 1. The improvements which are
obtained by multiple time stepping on the coarse meshes, are more pronounced for the
separated flow case as shown in Fig. 17c. Here, the interaction process between the shock
and the separating boundary layers is considerably accelerated, and converged values of lift
and drag are obtained more rapidly, see Table 1.
Before examining details of the flow solution, it may be appropriate to look at the
influence of the far-field distance on the lift and the drag of the wing. For this purpose, the
distance to the far field is varied by omitting grid planes at the far field. The effect of deter-
mining the free-stream quantities by superimposing the flow field of a horse-shoe vortex to the
onset flow is presented in Fig. 18. For far-field distances of 7 half spans, the results with and
without vortex differ only by 0.3%. The variation of the distance to smaller values shows,
that the use of the vortex reduces the sensitivity of lift and drag to the distance by about 50%.
In the results given below, the distance to the far field is kept at 7 half spans, and the horse-
shoe vortex terms are included.
The pressure distributions along several spanwise stations of the wing are displayed for
the attached-flow case in Fig. 19. The results of the 289x65~49mesh agree well with those
from the coarser mesh and with experimental data [20]. From Table 1, it is evident, that the
global features of the flow converge rapidly. Indeed, plots of the solution after 50 multigrid
cycles, which are displayed in Fig. 20, show virtually no differences to the fully converged
solution.
Results for the separated-flow case are displayed in Fig. 21. With the B.-L. turbulence
model, large discrepancies between computation and experiment occur. The size of the
separated region is underpredicted, and, consequently, the shock is located too far down-
stream. The agreement between computation and experimental data is greatly improved with
the J.-K.-mod. model, Fig. 22. The size of the separated region and the location of the shock
compare well with the data. The wall streamlines show a mushroom-type structure which is
typical for wings at high angle of attack. Table 1 shows, that approximately 80 multigrid
cycles on the fine mesh are sufficient for convergence to engineering accuracy. If the
193x49~33grid is considered to be sufficiently fine, 80 multigrid cycles require less than 40
minutes CPU-time on the Cray-2 computer.

Concluding Remarks
A cell-vertex scheme for the Navier-Stokes equations, which is based on central difference
approximations and Runge-Kutta time stepping, has been described. Using local time step-
ping, implicit residual smoothing with locally varying coefficients, a multigrid method, and
carefully controlled artificial dissipative terms, very good convergence rates are obtained for a
wide range of two- and three-dimensional flows. Details of the acceleration techniques, which
- 19 -

are important for convergence on meshes with high-aspect-ratio cells, have been discussed. In
general, engineering accuracy is obtained within 40 multigrid cycles on the fine mesh for
two-dimensional flows, and within 50-80 multigrid cycles in three dimensions. The accuracy
of the code is examined by grid refinement studies and comparison with experimental data.
For an accurate prediction of turbulent flows with strong separations, a modified version of
the nonequilibrium turbulence model of Johnson and King is introduced, which is well suited
for an implementation into three-dimensional Navier-Stokes codes. It is shown that the solu-
tions for three-dimensional flows with strong separations can be dramatically improved, when
a nonequilibrium model of turbulence is used.

Acknowledgement
The present work evolved within a scientific exchange program between DFVLR and NASA,
while the author stayed at NASA Langley Research Center. The author gratefully ack-
nowledges numerous discussions with M.D. Salas, Dr. R.C. Swanson, Dr. V.N. Vatsa. He is
grateful to Dr. R. Abid for sharing his knowledge about the Johnson-King turbulence model.
Furthermore, the author is indebted to C. Rossow for providing the cell-vertex Euler code as
the basis of the present work, and to B.W. Wedan for generating the three-dimensional grids.

Langley Research Center


National Aeronautics and Space Administration
Hampton, VA 23665
December 15, 1988

References
[l] Jameson, A.; Schmidt, W.; Turkel, E.: Numerical Solutions of the Euler Equations by
Finite Volume Methods Using Runge-Kutta Time- Stepping-Schemes. AIAA Paper 8 1-1259,
( 1981).
[2] Jameson, A.: Multigrid Algorithms for Compressible Flow Calculations. MAE Report
1743, Princeton University, Text of lecture given at 2nd European Conference on Multigrid
Methods, Cologne, Oct. 1985, (1985).
[3] Martinelli, L.: Calculations of Viscous Flows with a Multigrid Method. Ph.D. Disserta-
tion, MAE Department, Princeton University, (1987).
141 Radespiel, R.; Swanson, R.C.: An Investigation of Cell Centered and Cell Vertex Mul-
tigrid Schemes for the Navier-Stokes Equations. AIAA Paper 89-0548, (1989).
[5] Jayaram, M.; Jameson, A.: Multigrid Solution of the Navier-Stokes Equations for Flows
over Wings. AIAA Paper 88-0705, (1988).
[6] Johnson, D.A.; King, L.S.: A New Turbulence Closure Model for Boundary Layer Flows
with Strong Adverse Pressure Gradients and Separation. AIAA Paper 84-0175, (1984).
- 20 -

[7] Rossow, C.-C.: Berechnung von Stromungsfeldern durch L6sung der Euler-Gleichungen
mit einer neuen Finite-Volumen Diskretisierungsmethode. Dissertation TU-Braunschweig
1988.
181 Whitfield, D.L.; Janus, J.M.: Three-Dimensional Unsteady Euler Equations Solution Using
c
Flux Vector Splitting. AIAA Paper 84-1552, (1984).
[9] Klunker, E.B.: Contribution to Methods for Calculating the Flow about Thin Lifting
Wings at Transonic Speeds - Analytic Expressions for the Far Field. NASA TN D-6530,
(1971).
[ 101 Paisley, M.F.; Hall, M.G.: Improvements in the Formulation and Numerical Solution of
the Euler Problem for Swept Wings. Symposium Transonicum 111, DFVLR-AVA Gottingen,
May 24-27 1988.
[ l l ] Jameson, A.: The Evolution of Computational Methods in Aerodynamics. J. Appl.
Mech., Vol. 50, (1983).
[12] Baldwin, B.S.; Lomax, H.: Thin Layer Approximation and Algebraic Model for
Separated Turbulent Flows. AIAA Paper 78-257, (1978).
[13] Abid, R.; Vatsa, V.N.; Johnson, D.A.; Wedan, B.W.: Prediction of Separated Transonic
Wing Flows with a Non-Equilibrium Algebraic Model. AIAA Paper 89-0558, (1989).
[14] Abid, R.: Extension of the Johnson-King Turbulence Model to the 3-DFlows. AIAA
Paper 88-0223, (1988).
[15] Stock, H.W.; Haase, W.: The Determination of Turbulent Length Scales in Algebraic
Turbulence Models for Attached and Slightly Separated Flows Using Navier-Stokes Methods.
AIAA Paper 87-1302, (1987).
[16] Swanson, R.C.; Turkel, E.: Artificial Dissipation and Central Difference Schemes for the
Euler and Navier-Stokes Equations. AIAA Paper 87-1 107, (1987).
[17] Cook, P.H.; McDonald, M.A.; Firmin, M.C.P.: Aerofoil RAE 2822 - Pressure Distribu-
tions and Boundary Layer and Wake Measurements. AGARD-AR-138, (1979).
[181 Holst, T.L.: Viscous Transonic Airfoil Workshop Compendium of Results. AIAA Paper
87-1460, (1987).
[19] Harris, C.D.: Two-Dimensional Aerodynamic Characteristics of the NACA 0012 Airfoil
in the Langley 8-Foot Transonic Pressure Tunnel, NASA TM 81927, (1981).
[20] Schmitt, V.; Charpin, F.: Pressure Distributions on the ONERA-M6-Wing at Transonic
Mach Numbers. AGARD-AR- 138, (1979).
- 21 -

.,
c3
U
CL

V
I
I

I
-.J

I
m

I
A

I ;
H
Y

I ;
H
!
Y
I

4
- 22 -

Fig. 1 : Control volume around point (iJ,k). Fig. 2 : Control volume in a plane grid.

Fig. 3 : Definition of horse-shoe vortex Fig. 4 : Definition of variable d in


with respect to Cartesian coordinates. separated flow.
- 23 -

multigrid cycles C D=~ 0.0327


X,/C = 0.81

-1.5 -

-1 .o -
%
-.5 -
0-

-
.5-

1.0-

1.5, -.05L 1 1 1
0 .2 .4 .6 .8 1.0 0 .2 .4 .6 .8 1.0

Fig. 5 : Convergence histoyy, streamlines, and distributions of pressure and skin friction for
laminar flow around NACA 0012 airfoil, w4.5,do, Re,=5000, grid 257x65.
- 24 -

-2
M
0
*
-4

-6

-8

-lot \
a) Influence of fourth-difference dissipation on convergence, single time step on coarse
meshes, mesh 385x65, B.-L.model

2r c1 'r

-8

-lot
t
-12\
o zoo 400 600 aoo 1000 0 200 400 600 800 1000
mu1tigrid cycles multigrid cycles
b) Convergence for two time steps on c) Convergence for nonequilibrium J.-
coarse meshes, K(4)=1/32, mesh 321x65, &mod. turbulence model, two time steps
B.-L.model on coarse meshes, Id4)=l/32, mesh 321x65
Fig. 6 : Convergence behavior for transonic flow around RAE 2822 airfoil, K 4 . 7 3 , ~ ~ 2 . 7 9 '
Re,=6.5x106).
- 25 -

.86

.84

.82
cl
.80

.78
<
C

.0140

.0135

.0130/
F 0
,0065

.0060

A Cd .0055
cd P *O1,? fl

-v
f
.0120 - .0050 A

.0115 0 .0045 0

b) Pressure drag c) Friction drag

Fig. 7 : Influence of grid density and artificial viscosity on global force coefficients for flow
around RAE 2822 airfoil @&,,=0.73, cc=2.79', Re,=6.5x106), computed with B.-L.model.
- 26 -

I
I

0 7 4 6 8 1 0 0 2 4 6 8 1 0 0 2 4 6 8 1 0
x/c x/c x/c

I
12 _x ' 0-3 12 x lo-' 12 x ' 0-3

10- 10-

8- 8-
c,
6- 6-
I

4- 4-

7- 2-
I
0 0

I I I I I 1 2--
n 3 4 6 8 1 0 0 2 4 6 8 1 0 0 2 4 6 8 1 0
I x/c x/c x/c

I a) 193x33 grid b) 385x65 grid c) 577x97 grid

Fig. 8 : Distributions of pressure and skin friction for different grid densities (RAE 2822 air-
foil, k 4 . 7 3 , a=2.79', Re,=6.5x106, 1~(~)=1/64),computed with B.-L. model .
- 27 -
-1 .t - 1 51-

-1 .c
cP
-.5

.5

1 .o

1.5 -1l.5L 1 5 1
0 .2 .4 .6 .B 1.0 0 .2 .4 .6 .8 1.0
x/c X/C

- --
12 x 1 0 - 3 experiment [ 171
0-3
r C, = 0.803, CD = 0.0168

10-

8-
cf
6-

4-
2tI-
I
' !
!

,
-2i I I 1
-2 I I I
0 .2 C .fi .8 1.0
0 -2 .4 .6 .8 1.0
x/c x/c
a) B.-L.model b) J.-K.-mod. model
CL = 0.849, CD = 0.0175 CL = 0.814, CD = 0.0162

Fig. 9 : Effect of the turbulence model on distributions of pressure and skin friction for RAE
2822 airfoil (M,=0.73, a=2.79', Re,=6.5x106, grid 481x97).
- 28 -
-1.5 -1.5)-

-1 .o -1 .o
cP cP

-.5 -.5

0 0

.5 .5

1 .o 1 .o

1.5 - 11.5L0 .2 .4 .6 .8 1.0


0 .2 .4 .6 .8 1.0
x/c x/c

12 experiment [ 171
CL = 0.743,C D = 0.0242

10

8
Cf

-2 - - -2L 1
0 .2 .4 .6 .8 1.0 0 .2 .4 .6 .8 1.0

a) B.-L. model b) J.-K.-mod. model


CL = 0.85 1, CD = 0.0289 C, = 0.810, CD = 0.0264

Fig. 10 : Effect of the turbulence model on distributions of pressure and skin friction for RAE
2822 airfoil (M-3.75, a=2.8lo, Re,=6.5x106, grid 481x97).
- 29 -
-1.5 -1Sk
-1 .o -1 .o
- CP %
-.5 -.5

0 0

.5 .5

1 .o

1.5 1.5L ,
0 .2 .4 .6 .8 1.0
x/c

12 10-5 experiment [ 191 D-J


CL= 0.390,CD = 0.0331
1 1
10

8
Cf
6

0 0-

-2 - - -2.
0 .2 .4 .6 .6 1.0
x/c x/c
a) B.-L. model b) J.-K.-mod. model
CL = 0.466, CD = 0.0434 CL = 0.386, CD = 0.0378

Fig. 11 : Effect of the turbulence model on distributions of pressure and skin friction for
NACA 0012 airfoil ( K 4 . 7 9 9 , ot=2.26', Re,=9.0x106, grid 321x65).
- 30 -

0.8

0.6

Experiment W I ,
CL 0.4 corrected
- B.-L. model
--- J.-K.-mod. model
0.2

Fig. 12 : Comparison of lift coefficient versus angle of attack for the NACA 0012 airfoil,
M,=€').7, Re,=9.0x106, grid 321x65.

0.8

0.6

CL 0.4

0.2

0 COl 0.02 0.03 0.04 0.05 0.06 0.07

Fig. 13 : Comparison of lift versus drag polars for the NACA 0012 airfoil, x 4 . 7 ,
Re,=9.0x106, grid 321x65.
- 31 -
2

0
-e4
-c1
c eQ -2
-
m
-
M -4
d
0

-6

-8

-1 0

-1 2
o zoo 400 600 aoo 1000
multigrid cycles

*"ti
1 .o

2.0L P
0 .2 .4 .6 .8 1.0 0 .2 .4 .6 .8 1.0
x/c x/c

Fig. 14 : Convergence history and distributions of pressure and skin friction for turbulent
supersonic flow around NACA 0012 airfoil, k = 2 . 0 , ad'', Re,=10.0x106, grid 321x65, com-
puted with the B.-L.model.
- 32 -

.
.

Fig. 15 : Mach contours and 321x65 mesh for turbulent supersonic flow around NACA 0012
airfoil, Mm=2.0, do, Rem=10.0x106,computed with the B.-L.model.
- 33 -

a) Partial view of 97x25~17mesh w

b) Spanwise sections of 289x65~49mesh at the wing tip

Fig. 16 : Coordinate grid around ONERA-M6 wing.


0
-
N
-
c1
9 -2
Q
-
m
- -4
M - \
\
0
4

-6 - 1/32

-8 - -8

-10 - -1 0

-1 2 1 1 u -1 2
0 200 4OOv 400 600 800 0 200 400 600 800
multigrid cycles multigrid cycles

a) Influence of fourth-difference dissipation b) Convergence for two time steps on


on convergence, single time step on coarse coarse meshes, €d4)=1/32, a=3.06', B.-L.
meshes, a = 3 . 0 6 O , B.-L. model model.

- I 1.28-

two time steps on


coarse meshes

-6 .6 -

-8

-10

-12 0 L 200 400 400 600


.T
.2 -

ob,-&
0 200 0 200 400 600
multigrid cycles multigrid cycles
c) Convergence for nonequilibrium J.-K.-mod. turbulence model, K(4)=1/32, ~ ~ 6 . 0 6 ~ .

Fig. 17 : Convergence behavior for transonic flow around ONERA-M6 wing K=0.84,
Re,=l 1.0x106, mesh 289x65~49
0.55 - ,- Mesh 137 x 29 x 25

CL 0.50

I I
0.45
0.055
\
With

CD 0.050

0.45
0 5 10
Farfield distance in half spans

Fig. 18 : Lift and drag coefficients versus distance to the far-field boundary for transonic flow
over ONERA-M6 wing M,=0.84, a = 6.06', Re,=11.0x106.
- 36 -
_.

.
8

1.2

-1 2-
L -
f
0 2
x/c
6 a 1.0

-1.2 -
1
0
I
2
I
4
X/C
I
6
I
a
I
1.0

-0- - 2 8 9 ~ 6 5 ~ 4grid
9 -.8 -
--- 1 9 3 ~ 4 9 ~ 3grid
3
- 4 - -.4 -
9 0-
u" 0-

.4 - .4 -
.8 - .8 -

1.2 L 0I
L
2
u
6 8 1.0 4
1.2L I
0
I
2
I
4
I
6
I
a
I
1.0

wall streamlines isobars on upper wing s u r f a c m

Fig. 19 : Pressure distributions and wall streamlines for ONERA-M6 wing, %=0.84,
a = 3.06', Re,=11.0x106, Id4)=1/32, computed with the B.-L.model.
- 37 -

z/a = 0.44

-
T'
1.2L 1 I I I I I 1.2L
.8

I I I I I 1
0 .2 .4 .6 .8 1 .o 0 .2 .4 .6 .8 1 .o
x/c x/c

-1.2

-.E

-.4

4 G 0

z/a = 0.80
.4

.8

--
0 .2 .4 .6 .8 1 .o
x/c

Fig. 20 : Comparison of pressure distributions on ONERA-M6 wing, K 4 . 8 4 , a = 3.06O,


Re,=l 1.0x106, K(4)=1/32, B.-L.model, after 50 multigrid cycles with the fully converged
result, two time steps on the coarse meshes.
- 38 -
-1.2

-. 8

-_1

$ 0

.4

.8

1.2L I I I I I I 1.2 I I I I I I
0 2 4 6 8 1.0 0 2 4 6 8 1.0
x/c

1
-1 2 7
-1 2-
o 0 measurement [20]
-8-
- 2 8 9 ~ 6 5 ~ 4grid
9 -.a -
--- 193x49~33grid
--- -4-
-1-

0" u" 0-

-
0-

.4 - .4 -
I

.8 - .8 -

4x'A
1.2L
0 2
X/C

T wall streamlines isobars on upper wing surf=-

Fig. 21 : Pressure distributions and wall streamlines for ONERA-M6 wing, &=0.84,
a = 6.06O, Re,=11.0x106, K(4)=1/32, computed with the B.-L.model.
- 39 -
-1.2 -1.2

-.a -.8

-.4 '
-.4

u" 0 . $ 0

'
.4 '
.4

.a . .a .

1.2L I I I I I I
0 .2 .4 .6 .8 1.0 0 .2 .4 .6 .8 1.0

-1.2
-
-.a -
-.4 -
0-

.4 -

.8 -

0 .Z .4 .6 .8 1.0

Fig. 22 : Pressure distributions and wall streamlines for ONERA-M6 wing, M,=0.84,
a = 6.06O, Re,=l 1.0x106, K(4)=1/32,computed with the J.-K.-mod. model.
1. Report No. 2. Government Accession No. 3. Recipient’s Catalog No.

NASA TM- 101557


__
4. Title and Subtitle 5. Report Date

A C e l l - V e r t e x Mu1t i g r i d Method f o r t h e January 1989


Navier-Stokes Equations .. _. ____-__-__
6. Performing Organization Code

8. Performing Organization Report No.

10. Work Unit No.

9. Performing Organization Name and Address 505-60-01-01


1 1 . Contract or Grant No.
NASA Langley Research Center
Hampton, VA 23665-5225
13. Type of Report and Period Covered
12. Sponsoring Agency Name and Address
Technical Memorandum
N a t i o n a l Aeronautics and Space A d m i n i s t r a t i o n
1 4 . 5 n s o r i n g Agency Code
Washington, DC 20546-0001

__.____-__
15. Supplementary Notes

R . Radespiel: NASA-NRC Resident Research Associate a t Langley, now a t


t h e I n s t i t u t e f o r Design Aerodynamics , DFVLR, Gottingen, West Germany (FRG).

_-
is. Abstract
A c e l l - v e r t e x scheme f o r t h e Navier-Stokes equations, which i s based on c e n t r a l
d i f f e r e n c e approximations and Runge-Kutta t i m e stepping, i s described. Using l o c a
time stepping, i m p l i c i t r e s i d u a l smoothing, a m u l t i g r i d method, and c a r e f u l l y
c o n t r o l l e d a r t i f i c i a l d i s s i p a t i v e terms, v e r y good convergence r a t e s a r e o b t a i n e d
f o r a wide range of two- and three-dimensional f l o w s over a i r f o i l s and wings. The
accuracy of t h e code i s examined by g r i d refinement s t u d i e s and comparison w i t h
experimental d a t a . For a n a c c u r a t e p r e d i c t i o n o f t u r b u l e n t flows w i t h s t r o n g
separations, a m o d i f i e d v e r s i o n o f t h e n o n e q u i l i b r i u m t u r b u l e n c e model of Johnson
and King i s i n t r o d u c e d , which i s w e l l s u i t e d f o r an implementation i n t o t h r e e -
dimensional Navier-Stokes codes. I t i s shown t h a t t h e s o l u t i o n s f o r t h r e e -
dimensional f l o w s w i t h s t r o n g s e p a r a t i o n s can be d r a m a t i c a l l y improved, when a
n o n e q u i l i b r i u m model of t u r b u l e n c e i s used.

____ - - . .._ __ ___--__-


17 Key- Words (Suggested by Authoris)) 78 Distribution Statement
N a v i er-Sto kes Equations Unclassified-Unl i m i t e d
Numerical Methods S u b j e c t Category 02

- 22. Price
19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of pages

Unc 1ass if ied Unclassified 40 A0 3

View publication stats

You might also like