NASA Technical Reports Server (NTRS) 19930018276: A robust quasi-simultaneous interaction method for a full potential flow with a boundary layer with application to wing/body configurations

16 Pages·0.72 MB·English

by  NASA Technical Reports Server (NTRS)

A ROBUST QUASI-SIMULTANEOUS INTERACTION METHOD FOR
A FULL POTENTIAL FLOW WITH A BOUNDARY LAYER WITH
APPLICATION TO WING/BODY CONFIGURATIONS *)
A.J. van der Wees **) and J. van MulJden ***)
N68-274,6
National Aerospace Laboratory
P.O. Box 90502; 1006 BM Amsterdam
The Netherlands
Tel: [phone]; Fax: [phone]
Abstract formulated. Robustness of the interaction algorithm
has been formulated as the main requirement. Next,
The MATRICS flow solver calculates the the basic concepts of the vlscous-inviscld
inviscid transonic potential flow about a winE/body interaction solver will be given. Subsequently, a
semi-configuration. At present, work is in progress short system overview will be given of the MATRICS-
to extend MATRICS to take viscous effects into V flow solver, followed by an analysis of the
account through coupling with a boundary layer boundary layer equation system and on the inter-
solver. This solver, MATRICS-V, is based on robust action law. Finally, some preliminary computational
calculation methods for the boundary layer, the results will be presented.
outer wing flow and their interaction. MATRICS-V is
intended for (inverse) design purposes. The 2. Startln_ vosltion for the
boundary layer and wake are based on an integral development of MATRICS-V
formulation of the unsteady first order boundary
layer equations, the inviscid method is the The MATRICS-V vlscous-lnviscid interaction
existln 8 MATRICS potential flow solver and the solver is a follow-up to the MATRICS three-
interaction algorithm is of the quasl-slmultaneous dimensional Inviscid transonic potential flow
type. solver. This flow solver solves the full potential
equation in strong conservation form on a grid of
The paper gives a progress report on the C-H or C-O topology. This grid is generated using
coupled potentlal-flow boundary-layer method for the MATGRID grid generator (Tysell and Hedman,
transonic wlng/body configurations. Ref. 2). The solver uses a fu!ly conservative
finite volume dlscretlzatlon and a multlgrld
i. Introduction solution method. The dlscretlzatlon scheme is
second order accurate in the mesh size in subsonic
Computation methods for three-dlmenslonal parts of the flow, and first order accurate in
transonic potential flow are an important component supersonic parts of the flow. For the capture of
in design systems for civil aircraft wing/body supersonlc/subsonlc shock waves a Godunov-type
configurations. Accurate performance prediction shock operator is used. Options for fully-
under these conditions (e.g. lift-drag analysis) conservative as well as non-conservatlve shock-
requires that the transonic potential flow solver capture are available. Details on the MATRICS flow
is coupled with a boundary layer solver to account solver can be found in references 3, 4.
for the viscous effects on the wing pressure MATRICS provides data for a lift-drag-dlvlng-moment
distribution in a sufficiently accurate way. For analysis. Substantial research has been performed
some years NLR has available its own developed on the reliable prediction of drag by the MATRICS
system MATRICS (Multi-component Aircraft TRansonic flow solver. Findings of this research have been
Inviscld Computation System) for the calculation of reported in reference 5.
the three-dlmensional invlscld transonic potential
flow about wing/body configurations (Ref. I). The The development of MATRICS-V is partly based
MATRICS-derlvatlve MATRICS-V -now under develop- on experience at NLR with the modelling of two-
ment- will calculate the interaction of the dimensional strong viscous-invlscld interaction on
invlscld potential outer flow and a viscous airfoils (Refs. 6, 7).
boundary layer on the wing of a transonic transport
winK/body configuration. The ultimate objective of 3. Reoulrements
the development of MATRICS-V will be the embedding
of this system in a wing design system, also MATRIGS-V is designed to calculate the
currently under development at NLR. influence of the wing boundary layer and wake on
the lnviscid potential outer flow about a given
This paper gives a progress report on the transport-type wlng/body configuration from
development of MATRICS-V. Firstly, a description subsonic up to and including transonic cruise
will be given of the MATRICS three-diamnslonal conditions. Laminar as well as turbulent boundary
transonic potential flow solver, being the starting layers should be calculated, where the turbulent
point for the development of MATRICS-V. Secondly, boundary layer is allowed to be mildly separated.
the requirements for the development of the The MATRICS-V method should be robust; a guaranteed
viscous-inviscld interaction solver will be converged solution should be obtained for realistic
flow conditions. The method should also be about
ten times faster than a three-dlmenslonal Reynolds-
*) This research has been performed under averaged Navler-Stokes solver, to enable wing
contract with the Netherlands Agency for design applications.
Aerospace Programs (NIVR 01802N)
**) Research scientist, Num. Math. and Appl. Frog. &. Basic concepts
Dep., Informatics Division
***) Research Scientist, Theoretical Aerodynamics The Coupling of an inviscid outer flow with a
Dep., Fluid Dynamics Division viscous boundary layer flow will be done using an
integral method formulation of the boundary layer 5. Descriotion of the calculation metho d
equations. With integral methods there is no
explicit formulation of a turbulence model, but the 5.1 Inviscid method
system of equations is supplemented by suitable
empirical closure relations. The integral method is The full potential equation
reasonably easy to implement, because in the
integral method the three-dimensional flow problem (pub -0, (11
is formulated as a two-dimensional problem on the axi
wing surface and the wake center-surface. This is a
!
great advantage in code development.
ui = u®i + --a_ , (2)
Until a few years ago, the boundary layer ax i
equations have always been used in their steady
form to compute a steady boundary layer flow. The
main advantage of this formulation over an unsteady q2 = (uZ) 2 + (u2) 2 + (u3)2 , (3)
formulation has always been its lower computation
time. However, using the boundary layer equations
in their unsteady formulation, integrating them
towards a steady solution, the boundary layer
method is particularly well suited for
vectorization and hence for implementation on
todays supercomputers, see Van Dalsem and Steger is solved using a fin£te voium-e discretizatt6_ :_
(Ref. 8), Swafford and Whitfield (Ref. 9). An even formulated on a curvIilnearc6brdinate sys_em_and a
more important advantage of using the unsteady multigrid method employing ILU/$1P smoothlng. At
boundary layer equations for solving steady present only semi-conflgurations can be considered.
boundary layer flow is the robustness of this
approach. With the unsteady formulation a simple The boundary conditions arei
time-integration scheme will in any case produce an • on wetted surfaces u. - qS; (51
i
unsteady answer. With the steady formulatlon a for the inviscid flow solver, the source
space-marchlng scheme has to be used and the strength S - 0 on the body, els e S is computed
specification of initial data proves to be more by the boundary laye£solver; .....
difficult in a space-marchlng scheme than in a • in the sylaetry plane u, - O, (6)
time-integratlon scheme. Also in case the steady • in the far-field_eXcept downstream_ _=- 01:(7_ m
boundary layer solution is non-unlque or non- • in the far-field downstream (Trefftz-plane)
existent, the solution produced by an unsteady
formulation is probably more useful, than the a2_ = 0, (8)
solution produced by a steady formulation.
Two candidate interaction algorithms have been
considered, namely the seml-lnverse and quasi- where _I is the chordwise (wrap-around) grid
simultaneous algorithm. (The more sophisticated coordinate,
simultaneous algorithm is too expensive to • across the prescribed vortex sheet z
implement in MATRICS because of the fully implicit
relaxation algorithm employed in the existing [q] " q" - q-, (9)
inviscid potential outer flow solver). Experience bua] = (pqS)" - (pqS1- ;
by Ashill e.a. (Ref. i0), pp. 35, and Cebecl, Chen
e.a. (Ref. 111 indicates that the seml-inverse for the inviscid flow solver the Jump acros_
algorithm lacks robustness in difficult flow cases. the wake, (.)+ f (.)', Is computed by the
Therefore the quasi-simultaneous interaction boundary layer solver,
algorithm is used, being the best available • across the C-O topology branch cut that
interaction algorithm that can be implemented in extends from the tip section to the far-field
interaction with the inviscld potential outer flow lateral boundary
solution algorithm.
b] - o, b_] - o. (lO)
In the inviscid potential outer flow solver
the "blowing velocity approach" is used to account
for the boundary layer effects on the inviscid 5.2 Viscous method
outer flow, which means that a source strength is
specified at the body surface and a Jump in source The steady first order boundary layer
strength at the wake. This approach is reasonably equations, describing conservation of mass and
standard, while experience by Chen, Li e.a. momentum in a general rlght-handed coordinate
(Ref. 12) reports that this approach is more system, can be found in Myring (Ref. 15). Adding
reliable than the boundary layer displacement time-dependent terms and using first order integral
approach. thicknesses, the boundary layer equations can be
integrated with respect to z (normal to the wing
Experience by Chow (Ref. 13) indicates that it surface or wake centerline), using the momentum
is mandatory to prescribe the pressure Jump across equation in z-direction to eliminate the pressure.
the wake as a boundary condition to the outer Then the integral equations are obtained in the
inviscid flow. The latter Jump influences the Kutta form
condition prescription in Uhe outer inviscid flow.
and this is essential in order to compute realistic x-mom4nt-um:
lift values for the wing. The pressure Jump across
1 a6_
the boundary layer is computed as in Chow (Ref. 13)
and Lock and Williams (Ref. 14). _z at -
In the latter equation (13) the instationary
entrainment coefficient is an extension of the
(Ii) unsteady two-dimensional definition as used by
h, ax [h, Jaxthd'k J"
Houwink (Ref. 16).
Subsequently an expression is needed for the
calculation of the non-dimenslonal mass flux S
through the surface of the wing and wake,
f 1 representing the displacement effect of the
boundary layer on the invlscid outer flow. Assuming
J an inviscid flow between the stream surface z - 0
h2 Oy l h2 _ ay J ay[h,J and the displacement surface _*, which is a stream
surface for the inviscid flow, the following
expression can be obtained from (13) by setting
6 -0:
"I
,.__1 .- ,J'
1
* e22k2 - _C_I;
z
(14)
y-momentum: this is equivalent to S - w/q, yielding the usual
interpretation of S in steady inviscid flow.
In the used body-conforming non-orthogonal
jar coordinate system the y-axle is in the spanwlse
direction of the wing, while the x-axis is in
chordwlse direction wrapping around the wing and
the wake.
Next, a streamline coordinate system is adopted, in
which the variables (now denoted with tildes)
reduce to their familiar form, see Myrlng
lae, + 821L2 ,la ÷ --l-a-[ ] + 13}÷ (12) (Ref. 15). Transformation to and from the stream-
hIax [ hl _ 8X J ax line coordinate system is done whenever necessary.
Thus it is possible to derive equations using the
well-known integral parameters in streamline
coordinates and in the curvillnear system.
The equations (II) to (14) are thus written in the
basic variables _n, H, q and C, where the cross-
+ 022 + --i + 12 + flow parameter C is defined as
hz ay [ h2 _ ay Jay
c - slgu(_.)-_= (15)
Reduction of the number of unknowns to the four
basic variables is established by prescribing a set
of turbulent velocity profiles in the streamline
coordinate system, while the density thickness 6_
is eliminated using the Crocco relation,
1c . prescribing a parabolic distribution between
÷eull" _ _2'
velocity and temperature (Ref. 17). For the time
being no velocity profiles are used, but proven
entrainment : closure relations taken from accepted two-
dimensional methods for attached as well as
separated flow (Ref. 16, 18, 19) have been
;qJl at -q at (,x[h, j implemented in a first version of the code. In a
next version, more physical closure relations and
velocity profile families will be used.
Initial conditions for the turbulent boundary
layer calculation are generated by the BOLA-2D
[h, JJ ax h,ay
solver (Ref. 20), which calculates the laminar
quasi-two-dimensional flow in the leading edge
region of the wing.
(13) Boundary conditions are set at the wing root, where
derivatives in spanwise direction are assumed to be
where q is the velocity, J is the Jacobien of the zero, and the wing tip, where zero lateral deriva-
transformation from physical to computational tives _n local sweep direction are prescribed, see
space, and overbars denote boundary layer edge Cross (Ref. 21). Far downstream, a zero gradient
values. The density thickness 6_ is the integral of condition in chordwise direction is specified for
(P'a)/P over the boundary layer. all quantities.
The
system of equations (Ii) to (13) is solved specified on the wing. Further downstream, a quasi-
in combination with an interaction law (to be simultaneous interaction algorithm is used. In this
discussed in the next section). The complete set of formulation the inviscid flow calculation is done
equations can be written symbolically as in direct mode with a prescribed source strength S
on the wing and the wake and a prescribed velocity
_% + Au. + BuT +Du = f, (16) jump across the wake. The viscous calculation is
q done in quasi-simultaneous mode with a prescribed
Invlscid wall velocity, which has been corrected
for boundary layer curvature effects. Thus the
where u = [Sn, H, q, C]. boundary layer is computed effectively with a i
prescribed edge velocity instead of the invlscid
The system (16)+_.appears _o be hyp?rbo!ic in wall velocity. Cebeci, Clark e.a. (Ref. 28) have
practice. Discretizatlon is done according to t_e shown that such an approach avoids the initiation
directions of the characteristlcs in (x,t) and of undesirable pressure fluctuations in the |!
(y,t)-space, using a matrix-split procedure trailing-edge region at reasonably large angles of
(Ref. 22). Thus equation (16) is discretized as attack. The boundary layer computation computes a
new source strength S and velocity Jump Aq, which =
are used as the subsequent input for the
+A' .u•A- • •B-%u f (17) interactive calculations. =
q
_, Analysis of the system of eouatlons
In smooth parts of the flow second order accurate _,i Analysis of boundary laver eouation system
differencing will be obtained using a scheme as for
example in Spekreijse (Ref. 23). To analyze the properties of the boundary
The system of equations (17) is solved using the layer equations formulated in chapter 5 (Eqs. (1),
fully implicit backward Euler time-lntegration (2). (3). (18)) we assume the following simplifying
scheme proposed by Stager and Warming _(Ref. 24) and conditions:
Yee (Ref. 25). • orthonormal coordinate system, i.e. hl-hz-l,
g-O; kl-kz-k3-O; ll-lz-13-O;
5.3 Interaction law • outer stre-mline aligned with the x-axis, i.e.
&/_-l._/_-0;
In order to avoid a breakdown of the boundary • closure conditions as in Cousteix and
layer formulation in separated flow regions an Houdeville (Ref. 20), i.e.
extra equation is needed which modifies the _=
inviscid flow boundary layer edge velocity q.
Usually a highly linearlzed form of the inviscld
out.e.r....potential flow is taken, for ex.ample the t_o:
dimensional Hilbert-lntegral formulation as used by
_,,--cZ(_-_.), _--c_. (19)
Veldman (Ref. 26). An even more simplified form is
given by Williams (Ref. 27). In this paper the
latter form is slightlymp_d.ified, but still derived With H-_/_x, and H1-(_-_1)/_ n we obtain an
from the linearized potential equation. This wiil equation system _
be discussed in more detail in section 6.2. In its
simplest form the interaction law can be written as Eu_ + Au, + Bu_ - f, u = [In0n.H.Inq.lnC].
q
8q _ '_ S = 8J--r- _----_Sk. (18)
8s _x [SsJ _Ax where
Time-dependent terms are obtained from the
instationary form of S (equation (4)). paying
special attention to the limiting case for M * 0.
Two remarks are made: H - _(IP-I); r 1 - H(M2-X) + _2&r(_l+l) 0
I. Equation (18) is written along streamlines,
which implies that the streamline dlrectlona E - S + HI - &r(_l+l) H1 + I - 1_2(R+Hl-#r(l_l) 0 J ,(20)
are known from the invlscld flow solver and
Ir + _. H(I-_ 2) + _lqlr(R+l) H
are kept fixed during a viscous calculation. o+2o+I
2. Equation (18) is a law in correction form,
indicating that it will not influence the
converged solution. In this way the
interaction law can be shown to be essential H1 H1 H1(1-k2) 0
* - ,(21)
to avoid breakdown of the boundary layer
formulation, but once convergence is obtained _i -Ki(#r+l) Kd -KLMq_r(H+I) 0
it does not affect the final solution.
0 2+N2 l
_.4 Vfscous-invtscid interaction aleorithm
H-1 #r+l (1-_2) fl+(l_+l)_lq| r+_2 -2 H-1 ]
The leadin S edge part of the boundary layer
will be calculated in direct mode using the program S _r+_ (l-_2)n + (_P,l)_q| r H
;(22)
BOLA-2D (Ref. 20). This presupposes that the flow s-
0 0 0 0
does not separate in this part of the boundary
layer. The inviscid flow computation in this part
is done in direct mode with a source strength S C(H-1) C(Ir+l) C((R-I) (2 +_2)+(_,l)_qir) 2C(H- 1)
the parameter e_O at f_-O and all values for _ are positive for all H, so that
now all characteristics originate from the upwind
(23a) direction.
2
For the case a/aye0 we find the following
characteristic directions in the x,y-plane
- 2 * (_-i)_z (23b) (A-dy/dx) :
• No interaction law (Ki-O , Kd-l):
H = (H+I) (I÷6"=) - i, (23c) A - 0, (27a)
H_ = dHJdH. (23d) = H-l, (27b)
Setting K_-I and Ki-._/(_Ax) in the third row of
eqs. (20)..(22) reproduces eq. (18). H/A i + ((H-I)H/-HI(I+6r)+l+6r) _ +
Following Myring, reference 16, the
characteristics of an equation + 1 + _r " 0; (27c)
Eut + Au_ + f - O, u " [uI ... un] as in the x,t-plane a characteristic direction
changes from upwind to downwind at separation
can be obtained from (i.e. at H_-O).
With interaction law (Ki_), Kd-1) in case Ks-I,
det(E-AA) - 0, A = dt (241 _2-0, c-O:
dx
I - O, (28a)
We first consider the quasi t_o-dimensional flow
case (a/ay-0) where no interaction law is used,
i.e. Ki-O , K_-I. This way we find the following (-(H+I)IH_+(H+I)HI)I "_+ ((2H3iHZ-2H-2)H/ +
expressions for the characteristic directions A:
A = H, (25a) +(-2Hi+I)HI+Hi+H+I)A l +
=_t#,(_+l) ' (25b) +((HZ-I)H/+HZ+H+2)A • 1 " 0; (28b)
two positive and one negative value for _ are
- H[A 2 + ((H÷I)H_-Hx÷I-(H÷I)6ztI_) A ÷
found in case of realistic Hl(H)-functions.
- _f +H_+(_+l)6:H,_. o. (25c) The characteristics in the x,t-plane for the
more general case with an interaction law (Ki*O)
and without the additional settings Kl-1, _[_-0 and
For realistic Hl(H)-functions, viz. with HI-H1/H ¢-0 can only be analyzed numerically. We find the
the roots of the latter quadratic equation are following expressions for the characteristic
real. Consequently, the equation system is fully directions :
hyperbolic. If H_<O the values for _ are greater
than zero. At separation, that is at H_-O, in the - H, (29a)
minimum of the HI(H) curve, one eigenvalue I passes
through zero, which means that the characteristic
a3A3 + alli + all + so - O, (29b1
direction changes from upwind to downwind there.
Consequently, the equation system models the
corresponding physical behaviour. where a3 - det A and sO - det E.
A well-defined interaction law has the property
We continue the analysis for the quasi two- that a3_O for all _. This gives the following
dimensional case (8/8y-0), but now for the case relation between Ki and Kd:
where an interaction law is used to solve for lnq,
-H_
i.e. Kt,_O, K_-I. In case Kt-I, M_-O and e-O we find _> K,,,
the following expressions for A:
(HI,I[ -t!1) (1t+1) -,f_qlt[ (H÷I) +Htl[ (1-M l) -Hl_l, (H+I)
I - O, (26a1
(30)
if H_<O the value Ki-,0 may be chosen (direct mode).
(26b)
_ mH, If H_<O _he value Ki-O.OIK d is suitable for all _z
and all _, in the sense that numerical co_ucation
of _ indicates that then for 0<M2<1.7 all values
((H+I)2H/-(H+I)HI) Az + (-H/(2Hil)(H+I) + for _ are real and positive, while in some rare
cases complex values are found with a positive real
part. The boundary layer equation system (and
+H I(2H+31 - 11 A + (H+I)(HH_-HI) - 0. (26c) interaction law) have consequently favourable
properties for use in the transonic flow regime at
For realistic Hl(_)-functions, namely with (H+I)H_- cruise conditions. Because all characteristic
Hi<O and HH_-HI<O three values for A are positive, directions generally originate from upstream,
one is negative. For small e<O it can be shown that initial conditions for Oil, H, q and C have to be
specified upstream,while downstreamno
initial Parameters are now introduced for the direct and
conditions need be specified. inverse parts of eq. (35):
6.2 Provertles of interaction law
Ke8 [laq""vKi(S-qS1kk)'s8Sa-_k] (38)
Following Lock and Wiliiams (Ref. 14), we
consider a source strength S approximately normal
to the surface, a velocity q approximately parallel The parameters Ki and Kd can be used to perform
to the surface, and define a perturbation potential direct, inverse and quasl-slmultaneous computations
as as follows: _ =_
direct : Ki - O, Ke - i, (39a)
A--_q- _,. AS - _,. (31) inverse : Ki - I, Kd - 0, (39b)
q quasl-slmultaneous: Ki > O, Kd - I. (39c)
The outer potential flow can be described by the With quasl-simultaneous computations,
linearized percdrbacion potential equation inspection Of equation (3.8) shows that the
inviscid outer flow solution is modified locally at
(1-Mz) _,, + _ss = O, places where S differs much from the previous
solution Sk. In the first =few iterations between
inner and outer solver chls will lead to non-
with 4z = AS at Z = O. (32) physical velocity distributions in the boundary
layer, especially at the wake where= % velocity
In the context of the definition of an interaction difference will occur between upper and lower side
law we will interpret• Aq/q and AS as corrections for lifting cases, showing the deficiency of not
with respect to a starting solution, i.e. modelling the circulation in the interaction law.
The inviscld flow solver will therefore have to
_q .q-q_, As -s-sk, (33) compute thls circulation effect on Its own.
q q
Conslderlng the== two-dlme_nslonal steady for_m of
where qk is given by the preceedlng outer flow S and using the property of mass conservation along
computation and Sk by the preceedlng viscous flow streamlines in the outer fl0W, we will use in
calculation. A solution to (32) in subsonic flow in eq. (38): |
Fourier space is given by
S = 86" (40) =B
. Cei,.e-D,s, $ =_/_'Z_" (34) 8S
Algebraic manipulation of (31), (32) and (34), It can now be observed that with Kl>O and Ks>0 in
using the derivative of Aq/q with respect to s in "equation (38) 8q/ge has a positive correlation with
order to obtain non-imaginary quantities, then 86"/8s, which is a desired property in subsonic
leads c0-ala-w_in=Fourler_space of c_e form flOW.
[!aq_ i aqk I . u(s_sk). (35) 7. Prelimlnarv comvutatlonal results
tq as q= _j
The MATRICS invisctd flow solver has been
tested extensively and appears to be a robust
In physical spats we may consider (35) as the method (Ref. 3, 4). In thls section attention will
leading term in an integral form of an interaction be given to the boundary layer solver wlth
law and use it in thls simplified form. interaction law and its robustness. Finally a fully
In supersonic flow only the rlght-running wave is converged solution of the whole vlscous-invlscld
considered, interaction computation will be shown.
Starting wlth an invlscld velocity distribution
= Ce1_,e_lhz. (36) provided by the MATRICS outer solver, some
calculations have been made to obtain a converged
Manipulation now yields • non-imaginary interaction boundary layer solution only.
law of the form Initial values for the boundary layer parameters
were obtained from a _wo-dlmenslonal flat plate
# q_qk = _ (S_Sk) " (37) solution. These calculations, without any inter-
q action with the invisctd outer solver, show _'le
changes made by the interaction law.
Considering eq. (37) as the leading term of an
integral form of an interaction in physical space, The first case is a NACA-O012 straight wing at
we now have a form of interaction law that couples a Mach mmber of 0.70, zero angle of attack and
S to q instead of 8q/gs, which cannot be Reynolds number of 9 million. Only the root section
incorporated in the system of boundary layer of the wing will be shown, because Chls section is
equations in • simple way. The subsonic law (35) • symmetry plane wlth a t_o-dlmensional flow by
appears however co produce useful results both in definition, and therefore allows comparison with
subsonic and supersonic flow, which is mainly due two-dimensional methods. In figure la the invlscld
Co the fact that strong interaction occurs in velocity distribution at the root section is given.
subsonic parts of the flow (at supersonic-subsonic It can be shown that the trailing edge stagnation
shock waves end at trailing edges). This relaxes behaviour causes • severe boundary layer growth,
the need for an explicit supersonic law, also leading to separation if no interaction law is used
considering the findings described in section 5.1. to adapt the velocity distribution. Figures lb, c,
We will therefore use eq. (35) both in supersonic d show the resulting velocity, momentum r_hickness
and subsonic flow. and shape factor H at Reynolds number of 9 million,
calculated with interaction factor K±-0.051- The history is shown in figure 5a. Finally, in figure
resulting velocity distribution shows smaller 5b the convergence history of the residuals in the
decelerating velocity gradients and a less severe MATRICS outer flow solver are given as a function
stagnation behaviour at the trailing edge. In order of the number of smoothings. A reduction in both
to obtain a physically relevant viscous solution, maximum and mean value residual of 2.5 orders of
the outflow, computed from this very first boundary magnitude has been obtained, which is promising for
layer solution, can be added to the full invlscid further development of the interaction algorithm.
flow solver for an adjustment of the inviscid flow
to the boundary layer effect. A fully converged 8. References
viscous-inviscid solution requires a number of such
interations between the inviscid flow solver and i. Vooren, J. van der, Wees, A.J. van der,
the boundary layer/interaction law computation. Meelker, J.H., MATRICS, Transonic Potential
Flow Calculations About Transport Aircraft.
The second case is the same straight wing at In: AGARD Conf. Proc. No. 412, 1986
Math number of 0.799, .angle of attack of 2.26 2. Tysell, L.G., Hedman, S.G., Towards a General
degrees, and a Reynolds number of 9 million. The Three-Dimensional Grid Generation System.
Inviscid velocity distribution at the root section ICAS-88-4.7.4, 1988
(Fig. 2a) will certainly cause separation on the 3. Wees, A.J. van der, A Nonlinear Multigrid
upper side of the wing. Starting from a very simple Method for Three-Dimensional Transonic
flat plate momentum thickness distribution, large Potential Flow. Ph.D. thesis, University of
changes must be expected in the velocity gradients. Technology Delft, The Netherlands, 1988
In this case, the resulting momemtum thickness at 4. Wees, A.J. van der, Impact of Multigrid
the trailing edge is much below the viscous- Smoothing Analysis on Three-Dimenslonal
inviscid converged value, due to the changes in Potential Flow Calculations. In: Mandel, J.,
velocity distribution, while the shape factor H is e.a. (eds.), Proc. 4th Copper Mountain Conf.
only showing separation Just behind the shock wave on Multigrid Methods, SIAM Proceedings, ISBN
(Fig. 2b, c, d). In the viscous-inviscid converged 0-89871-248-3, 1989
solution separation will probably occur from shock 5. Vooren, J. van der, Wees, A.J. van der,
to trailing edge (Ref. 29). The non-physical Inviscld Drag Prediction for Transonic
velocity jump across the wake (Fig. 2b) is due to Transport Wings Using a Full-Potential Method,
the asymmetry of the computed flow and the AIAA-90-0576, January 1990 (accepted for
simplicity of the interaction law. In the iterative publication in AIAA Journal)
process between the full inviscid flow solver and 6. Veldman, A.E.P., Lindhout, J.P.F., Boer,
the boundary layer/interaction law, however, the E. de, Somers, M.A.M., VISTRAFS, a Simulation
errors due to the simplicity of the interaction law Method for Strongly-Interacting Viscous
should disappear, resulting in a converged solution Transonic Flow. NIR MP 88061 U, 1988
and an inactive interaction law. 7. Houwink, R., Veldman, A.E.P., Steady and
Unsteady Separated Flow Computations for
Figures 3a and 3b show the convergence Transonic Airfoils. NLRMP 84028 U, 1984
histories of the residuals of the x-momentum 8. Dalsem, William R. van, Steger, Joseph L.,
equation, entrainment equation and interaction law Efficient Simulation of Separate_ Three-
of the boundary layer equation system for the Dimensional Viscous Flows using the Boundary-
testcases presented in figures 1 and 2. Using the Layer Equations, AIAA Journal Vol. 25, No. 3,
same computational parameters, both convergence 1987
histories show a robust convergence, which means 9. Swafford, T.W., Whitfleld, D.L., Time-
that a converged solution can be computed even when Dependent Solution of Three-Dimensional
the inviscid flow velocity distribution and the Compressible Turbulent Integral Boundary-Layer
starting solution for the boundary layer flow do Equations, AIAA Journal Vol. 23, No. 7, July
not at all fit together. 1985
i0. Ashill, P.R., Wood, R.F., Weeks, D.J., An
At the time of writing of the paper, work on Improved, Semi-lnverse Version of the Viscous
the interaction algorithm was making good progress. Garabedian and Korn Method (VGK), RAE Techn.
A first fully converged interacting solution has Report 87002, 1987
been obtained for the attached flow around the II. Cebeci, T., chen, L.T., Chang, K.C., Peavey,
NACA-O012 wing at a Math number of 0.70, angle of C.C., An Interactive Scheme for Three-
attack of 1.49 degrees, and a Reynolds number of 9 Dimensional Transonic Flows, Proc. 3rd Symp.
milllon. The boundary layer is computed with an on Num. and Phys. Aspects of Aerodyn. Flows,
increasing part of the invlscid flow velocity in pp. 11-39 - 11-49, 1985
the first few iterations, while underrelaxatlon is 12. Chen, L.T., Li, S, then, H., Calculation of
applied to the source strength tha_ is used as Transonic Airfoil Flows by Interaction of
input to the outer solver. No boundary layer Euler and Boundary-Layer Equations, AIAA-87-
curvature effects have yet been accounted for. 0521, 1987
Using this procedure the inner and outer flow are 13. Chow, R., Solution of Viscous Transonic Flow
smoothly adapted to each other. over Wings, Computers & Fluids Vol. 13, No. 3,
In figure 4 the initial inviscid velocity pp. 285-317, 1985
distribution (viz. without a boundary layer effect) 14. Lock, R.C., Williams, B.R., Vlscous-lnviscld
is shown, together with the final viscous velocity Interactions in External Aerodynamics, Prog.
distribution and the corresponding boundary layer Aerospace Sci. Vol. 24, pp. 51-171, 1987
variables. The difference between the invlscid and 15. Myring, D.F., An Integral Prediction Method
viscous velocity is small, except at the trailing for Threedimensional Turbulent Boundary Layers
edge, where the boundary layer displacement effect in Incompressible Flow, RAE Techn. Report
is appreciable. This testcase has been calculated 70147, 1970
with thirteen viscous-inviscid iterations in order 16. Houwlnk, R., Veldman, A.E.P., Steady and
to obtain a converged lift value. The convergence Unsteady Separated Flow Computations for
TransonicAirfoils, AIAA-84-161189,84
Second-order Discretizations of Hyperbolic
17.
Cebec£, T., Smith, A.M.O., Analysis of Conservation Laws, Math. Comp. Vol. 49, No.
Turbulent Boundary Layers. Academic Press, 179, pp. 135-155, July 1987
1974 24. Steger, J.L., Warming, R.F., Flux Vector
18. Green, J.E., Weeks, D.C., Brooman, J.W.F., Splitting of the Invlscid Gasdynamlc Equations
Prediction of Turbulent Boundary Layers and With Application to Finite Difference Methods,
Wakes in Compressible Flow by a Lag- NASA-TM-78605, 1979
Entrainment Method, ARC R&M 3791, 1973 25. Yee, H.C., Upwind and Symmetric Shock-
19. Coustelx, J., Houdeville, R., Singularities in Capturing Schemes, NASA-TM-89464, 1987
Three-Dimensional Turbulent Boundary-Layer 26. Veldman, A.E.P., A new Quasl-Simultaneous
Calculations and Separation Phenomena, AIAA Method to Calculate Interacting Boundary
Journal, Vol. 19, No 8, pp. 976-985, 1981 Layers, AIAA J, Vol. 19, pp. 79'85, 1981
20. Bruin, A.C. de, Boer, E. de, Users Guide of 27. Williams, B.R., Coupling Procedures for
BOLA-2D, the NLR-metbod for the Calculation of Viscous-lnvlscld Interaction in External
(Quasi) Two-dlmenslonal Boundary Layers, NLR Aerodynamics. In: Proc. 4th Symp. on Num. and
memorandum AI-82-014 U, 1982 Phys. Aspects of Aerodyn. Flows, Long Beach,
21. Cross, A.G.T., Calculation of Compressible CA, January 1989
Three-Dimenslonal Turbulent Boundary layers 28. Cebecl, T., Clark, R.W., Chang, K.C., Halsey,
with Particular Reference to Wings and Bodies, N.D._ Lee, K., Airfoils with Separation and
YAD Note 3379, November 1979 the Resulting Wakes, Proc. 3rd Symp. on Num.
22. Hartwich, P,M., Split Coefficient Matrix (SCM) and Phys. Aspects of Aerodyn. Flows, pp.
Method with Floating Shock Fitting for 2-13 - 2-26, 1985
Transonic Airfoils. Proc. Int. Conf. on Num. 29. Holst, T.L., Viscous Transonic Airfoil
Meth. in Fluid dynamics, Oxford, July 1990 Workshop Compendium of Results. AIAA-87-1460,
23. SpekrelJse, S., Multlgrid Solution of Monotone 1987
i
_-=
__ffi
VtSU3NDLR ,.o4 1H1,H11,,1,I1,1] ,,H,,,/,lH/,,, ,,1,H1H,,1H1H1,,,1,,1,,1,,1,,1,,1,,,H,1,H11,,1H,,,,,1,,1,,,,-,'1,11111111
EQ.[f_.FLOW VI_|T'I' AT
1,60_
1.33.] 1.:rid
1._.2
0.00.1 o.ao.l
o._3_
o._,1
I
o.oo] 0,_.1
o'.0o...... o'.,,0 .... o'.eo .... _'._ " " " I'._ " " "2'.oo o'_ ...... o'._ ..... o'._ ...... I"._ ....... _;lO .... 2".o¢
X-AXIS X-AX%$
a) Inviscid flow velocity distribution b) Boundary layer tnvtscld edge velocity
NLR ,//Hill s//s/H, H/l/s/l/l,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
VISU3D _.o4 ,,,H,H ,,,,,,,, ,,,,,/H,H///,,,,H//////////////,,////
Xt0"2 STRENIV.MOM,DXSPL.THXCX
t.00_
0._3j
3.oo_
0.5o_ 2.so_
21
1,S0 _L ,,.
0.o0: 1.00J
o'.0o .... o'._ "" o'.im..... t'.2o.... £._ .... ':,'.o(
X-AXXI X-AXII
c) Boundary layer streumise Iolaentua thickness d) Boundary layer shape factor
Fig. 1 Boundary layer results for root section of a NACA-0012 straight wing at _-0.70. a=0, Re=9*106
NLR NLR
VISU3D ,.o4 VISU3( 1.04 J/l///I//Jl/ilt"_/ I//t#i/ I/I/I///J///I/I///l//iJ////l/i/l///I///J////I/J//l/l///J/l/l,l/l/l/J//i/l/lIllllllllll
EO.INV,FLOW VELOCITY AT INVIEIO li.L.8_l VlLgC?
I.G( I.W.
1.3."
i
0.8(
0.3_ o._z]
o.zT.]
0.,On,.1
oLOO o..,_ o.eo _.2o O.'W O;W ...... t;_ _'._m _".00
W-AXI|
a) Inviscid flow velocity diatrtbution b) Boundary layer inviscid edge velocity
//I/////////////////s//_//////////////////,
NLR I///'/"/:_/''/j
VISU3O I.o__..... "_'_.L:;" ,sI//_////////////////////////////////////#
|
4.00._
i
0.03 3._0_
J
0.67 3._.I
]
0.30
0._ .00.]
0.17
.
0.00 1,0_]
o'od .... o._o o.a_ 1.so _.so 2"0o
X-/_I$ X-AXIS
c) Boundary layer streamwile mom_ntu_ thickness d) Boundary layer shape factor
F_g. 2 Boundary layer_ results for root section of a NACA-0012 straight wing at M_ = 0.799, _ = 2.26,
Re = 9 _ 106