Table Of Content16th Computational Fluid Dynamics Conference, 23{26 June 2003, Orlando, FL
Failure of Anisotropic Unstructured Mesh
Adaption Based on Multidimensional Residual
Minimization
(cid:3) (cid:3)
William A. Wood and William L. Kleb
NASA Langley Research Center, Hampton, VA 23681
Anautomatedanisotropicunstructuredmeshadaptationstrategyisproposed,imple-
mented,andassessedforthediscretizationofviscous(cid:13)ows. Theadaptioncriteriaisbased
upon the minimizationof the residual(cid:13)uctuations of a multidimensionalupwind viscous
(cid:13)owsolver. Forscalaradvection,thisadaptionstrategyhasbeenshowntousefewergrid
points than gradient based adaption, naturally aligning mesh edges with discontinuities
and characteristic lines [Wood, AIAA 99-3254]. The adaption utilizes a compact stencil
and is local in scope, with four fundamental operations: point insertion, point deletion,
edge swapping, and nodal displacement. Evaluation of the solution-adaptive strategy is
performed for a two-dimensionalblunt body laminar wind tunnel case at Mach 10. The
resultsdemonstratethat the strategysu(cid:11)ers froma lackofrobustness, particularlywith
regard to alignment of the bow shock in the vicinity of the stagnation streamline. In
general, constraining the adaption to such a degree as to maintain robustness results in
negligibleimprovementto the solution. Because the present method failsto consistently
or signi(cid:12)cantly improve the (cid:13)ow solution, it is rejected in favor of simple uniform mesh
re(cid:12)nement.
Nomenclature Superscripts:
T
A~; ~ Flux Jacobian in conserved or auxiliary variables ( ) Transpose
a A Speed of sound ( )v Viscous component
x;y
F Flux of conserved quantities ( ) Spatial component of a vector
I Identity matrix
^{;|^ Cartesian unit vectors Subscripts:
‘ Edge length ( )i Node number
M Upwinding matrix ( )T Triangle
n^ Outward-normalunit vector
Over-bars indicate linear averages and tildes indi-
P Pressure
cate average values assuming linear variation of the
S Area
2 parameter vector.
s Entropy: ds=d(cid:26) dP/a
(cid:0) Variables lacking explicit dimensions are in nondi-
T Temperature
mensional form, normalized by the free stream values
U Conserved variables
of density, velocity, viscosity, and temperature along
u;v Velocity components
V~ Velocity vector with a characteristic length, as appropriate.
W Auxiliary variables: dW =(ds;(cid:26)~dV~;dP)
Introduction
x;y Cartesian coordinates
1
HOMASetal. identifyviscous(cid:13)owunstructured
Z Parametervector: p(cid:26)(1;u;v;H)
Tmeshadaptationasacriticalenablingtechnology
(cid:11);(cid:12) Generalized wave speeds in (cid:24);(cid:17) directions
for the inclusion of high (cid:12)delity computational (cid:13)uid
(cid:13) Ratio of speci(cid:12)c heats
dynamicsinthevehicledesigncycle. Thedevelopment
(cid:24) Curvilinear coordinate on triangle edge 12
ofrobustautomatedunstructuredgridgenerationand
(cid:17) Curvilinear coordinate on triangle edge 23
adaption methods will directly reduce the time re-
(cid:7) Minimization functional
(cid:30);(cid:30)(cid:21) Fluctuation in conserved or auxiliary variables quired for both the initial grid generationfor complex
vehiclesandthesolutioncomputationforcomplex(cid:13)ow
(cid:26) Density
(cid:12)elds. Notable achievements in this (cid:12)eld of unstruc-
(cid:4) Functional weighting matrix
turedmeshadaptionhaverecentlybeendemonstrated
(cid:10) Integration element 2
byVenditti andDarmofal, usingadjoints tooptimize
(cid:3) 3,4
AerospaceEngineer,AerothermodynamicsBranch,Aerody- the mesh distribution, and by Habashi et al., using
namics,Aerothermodynamics,&AcousticsCompetency. AIAA
anisotropic feature adaption.
member.
Concurrent work in (cid:13)ow solvers for hypersonic
ThismaterialisaworkoftheU.S.Governmentandisnotsub-
jecttocopyrightprotectionintheUnitedStates. aerothermodynamic applications has shown promise
1of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
3 n^1
for nonlinear multidimensional upwind (cid:13)uctuation
splitting schemes as being more accurate and less
grid sensitive than traditional locally one-dimensional
5{7 n^2 ‘1
(cid:13)ux splitting schemes. Unstructured mesh adap-
tion based upon the minimization of the (cid:13)uctuation (cid:17) (cid:24)
splittingdistribution has shownsomeremarkablenat- 2
8,9 10
uralalignmentpropertiesforlinear andnonlinear ‘2
‘3
model problems, but has not been demonstrated for
the full system of equations for viscous (cid:13)uid dynam-
ics. Alignmentofthemesh tothe bowshockhasbeen
1 n^3
shown to be a critical factor in the accuracy of hy-
11
personic solutions, and the successful extension of Fig.1 Elementaltriangulardomainfor(cid:13)uctuation
thenaturalalignmenttendencyinherentin(cid:13)uctuation splitting.
minimizationtoahighspeed(cid:13)uiddynamicscapability
could be the desired breakthrough to the cost barrier
The viscous(cid:13)uctuation is assembledfrom contribu-
ofincludingcomputational(cid:13)uiddynamicsinthevehi-
tions associated with each node of the triangle,
cle design cycle.
3
strTuhctisurpeadpemrepshresaednatsptainonauptroomceastsedbaasneidsoturpopoinc uthne- (cid:30)v = (cid:30)vi and (cid:30)vi = ‘2iF~~v(cid:1)n^i (3)
i=1
minimization of distributions for a nonlinear multidi- X
mensional upwind (cid:13)uctuation splitting scheme. The The edge of the triangleopposite to node i has length
adaption is local in scope and thus convenient to par- ‘i and outward unit normal n^i.
allelize. The demonstration case is for perfect gas Expressingtheinviscid(cid:13)uxintermsofitsJacobian,
two-dimensional Mach-10 (cid:13)ow over a sting-mounted dF~ =A~dU, the inviscid (cid:13)uctuation can be decom-
entry capsule pro(cid:12)le. posed along two edges of a triangle, denoted (cid:24) for 12
and (cid:17) for 23, as
Flow Solver (cid:24) (cid:17)
(cid:30) =(cid:30) +(cid:30) (4)
The (cid:13)ow solver is a nonlinear monotone second-
order accurate multidimensional upwind (cid:13)uctuation (cid:30)(cid:24) = 1‘1n^1 A~~(cid:1)(cid:24)U; (cid:30)(cid:17) = 1‘3n^3 A~~(cid:1)(cid:17)U (5)
12 (cid:0)2 (cid:1) 2 (cid:1)
splittingschemeasintroducedbySidilkover withex-
tensions to viscous (cid:13)ows.13 The solver produces exact Atransformationtogauxiliaryvariables,dU =gUW dW,
convective solutions along characteristics on unstruc- simpli(cid:12)es the expression for the Jacobian,
tured gridswhere one edgeof eachce1l0l is alignedwith (cid:30) =U~W(cid:30)(cid:21) (6)
thecell-averagedtransportdirection. Thischaracter-
istic alignment property of the (cid:13)ow solver is the basis (cid:30)(cid:21) =(cid:30)(cid:21)(cid:24)+(cid:30)(cid:21)(cid:17) = (cid:11)(cid:1)(cid:24)W (cid:12)(cid:1)(cid:17)W (7)
(cid:0) (cid:0)
fordrivingthemeshadaptionstrategydescribedinthe
with
following section. Traditional locally one-dimensional g g
1
(cid:13)ux splitting solvers do not share this method-of- 1 0 0 a2
1
characteristics property. u 1 0 a2u
While a full description of the second-order accu- UW =2 v2 0 1 a12v 3 (8)
V 1 T0
rateviscous(cid:13)owsolverisavailableinRef.6,including 6 2 u v (cid:13)(cid:0)1 T 7
6 7
vcreerti(cid:12)izcaattiioonniasnbdrivea(cid:13)lyidsakteiotcnherdesuhletrse, ttohep(cid:12)rorvstid-oerdcoerntdeixst- (cid:11) = 1‘1n^41 ~~; (cid:12) = 1‘3n^53 ~~ (9)
2 (cid:1)A (cid:0)2 (cid:1)A
for the residual derivatives that follow in the adapta-
and
tion section. The equations of motion are written in
vector form as 0 0 0 0
0 0 0 ^{
Ut = ~ F~ +~ F~v (1) A~ =V~I +20 0 0 |^3; (10)
(cid:0)r(cid:1) r(cid:1) 2 2
60 a ^{ a |^ 07
where F~ and F~v are the inviscid and viscous (cid:13)uxes 64 75
Upwinding is achieved through the introduction of
and U is the vector of conserved variables. Integrat-
an arti(cid:12)cial dissipation (cid:13)uctuation,
ingoveratriangularelement,Figure1,theright-hand
side (RHS) can be expressedvin terms of inviscid and (cid:30)(cid:21)0(cid:24) =sign(n^1 ~~)(cid:30)(cid:21)(cid:24); (cid:30)(cid:21)0(cid:17) =sign( n^3 ~~)(cid:30)(cid:21)(cid:17) (11)
viscous (cid:13)uctuations, (cid:30) and (cid:30) , as (cid:1)A (cid:0) (cid:1)A
which is linearly combined with the basic (cid:13)uctuation
RHS = ~ F~ d(cid:10)+ ~ F~vd(cid:10)=(cid:30)+(cid:30)v (2) to complete the inviscid contributionto the nodal up-
(cid:0) (cid:10)r(cid:1) (cid:10)r(cid:1) dates.
Z Z
2of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
Adaption Strategy of (cid:7) are formed using the chain rule as
T
The present adaption strategy for two-dimensional
unstructured meshes with triangular elements is local @(cid:7)i = 1 @(cid:30)(cid:21)T (cid:4)T(cid:30)(cid:21)T+(cid:30)(cid:21)TT@(cid:4)T(cid:30)(cid:21)T
in scope and allows anisotropic stretching. Adaption @xi 2 T 2 @xi ! @xi
X (13)
cycles are performed as sequential sweeps over the 4
domain performing the basic operations: node dele- +(cid:30)(cid:21)TT(cid:4)T@(cid:30)(cid:21)T ;
tion, edge swapping, nodal displacement, and node @xi3
insertion. The process is automated and a variety of
5
stopping criteria can be considered, such as a maxi- and
mum grid size, total number of cycles, orminimum to T
maximum residual ratio. @(cid:7)i = 1 @(cid:30)(cid:21)T (cid:4)T(cid:30)(cid:21)T+(cid:30)(cid:21)TT@(cid:4)T(cid:30)(cid:21)T
A node is marked for deletion if the (cid:13)uctuations, @yi 2 T 2 @yi ! @yi
a X (14)
bothinviscidandviscous, inall surroundingcellsare 4
belowathreshold. Thethresholdiscustomizable,and +(cid:30)(cid:21)TT(cid:4)T@(cid:30)(cid:21)T :
atypical choicewould be 2{3ordersof magnitude be- @yi 3
low the average nodal value. A trial reconnection of
5
the mesh is proposed without the node and the ele- Having de(cid:12)ned the gradient of the functional, the
ment (cid:13)uctuations are recomputed, from the existing method of steepest descent can be applied directly to
solution, to verify that the (cid:13)uctuations remain small. drive the nodal displacements.
If so, the node deletion is (cid:12)nalized. The weighting factor (cid:4)T is a symmetric positive-
de(cid:12)nitematrix. Weightingeachequationequallywith-
Anodeisaddedatthemidpointofanedge,creating
out regard to cell sizes results in (cid:4)T = I, while an
twomorecellsandthreemoreedges,ifthe(cid:13)uctuations 1
inthecellstoeithersideoftheedgeexceedathreshold. inverse area weighting yields (cid:4)T = STI. The deriva-
tives of (cid:4)T for these and other area-weighted choices
When comparingcell (cid:13)uctuations, the L2-norm of the
are available in Ref. 14. The present study is con-
(cid:13)uctuations from each of the governing equations are
ducted using equal weighting.
used,weightedbythe inverseofthe squarerootofthe
The discretecell (cid:13)uctuation, Eqn.7, can be equiva-
cell area.
lently written as
Edgesare(cid:13)aggedtobeswappedwhentherootmean
3
square(RMS) of the (cid:13)uctuations in the cells to either (cid:30)(cid:21) = 1 ‘jn^j ~~Wj; (15)
side are excessive. If the swapped edge maintains a 2 (cid:1)A
j=1
valid grid then the RMS of the (cid:13)uctuations is recom- X
puted,andtheswappedcon(cid:12)gurationisretainedifthe where Wj =WZ(cid:22)Zj. The summation is over the three
RMS has decreased. nodes de(cid:12)ning the triangle.
Nodal displacements are driven by the minimiza- The derivatives of the (cid:13)uctuation with respect to
tion of a functional that can be expressed in terms of themovementofonlyonenode, forexamplenode2of
either the conserved or auxiliary (cid:13)uctuations, as the the triangle, are
minimization of one implies the minimization of the
other. Because the goal of the nodal displacements is @(cid:30)(cid:21) = 1 ~y(W1 W3)
toachievealignmentofthemeshwithinvisciddisconti- @x2 22A (cid:0)
nuitiesandcharacteristiclines,theviscous(cid:13)uctuation (16)
is omitted from the displacement functional. +4 3 ‘jn^j @A~~ Wj +~~@Wj ;
Thenodaldisplacementfunctionalisconstructedfor j=1 @x2 A @x2 !3
each node i as X
5
and
(cid:7)i = 21 (cid:30)(cid:21)TT(cid:4)T(cid:30)(cid:21)T; (12) @(cid:30)(cid:21) = 1 ~x(W3 W1)
XT @y2 22A (cid:0)
(17)
wherethesummationisoveralltrianglesconnectedat +4 3 ‘jn^j @A~~ Wj +~~@Wj :
node i and (cid:4) is the weighting matrix. The derivatives j=1 @y2 A @y2 !3
X
5
The derivatives of the (cid:13)uctuation with respect to
aThesumoftheabsolute valuesoftheviscousdistributions movements of node 1 or 3 of the triangle follow di-
isused,asthenetcellularviscous(cid:13)uctuationisidenticallyzero. rectly as cyclic permutations.
3of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
The (cid:13)ux Jacobian in auxiliary variables, Eqn. 10, + + @Z2
(cid:11) (cid:12) WZ(cid:22)
has the approximate derivative T (cid:0) @x2
X(cid:0) (cid:1) + (27)
@ ~ @V~ = @(cid:11) WZ(cid:22)(Z1 Z2)
A I; (18) @x2 (cid:0)
@x ’ @x T (cid:20)
X +
@(cid:12)
and similarly for the y derivatives. It is assumed that WZ(cid:22)(Z3 Z2) ;
moving a node does not change the solution at the (cid:0) @x2 (cid:0) (cid:21)
other two nodes of the triangle, so that the variation and
of the cell-average (cid:13)ux Jacobian scales like one-third
(cid:0)1
the variation of the velocity at the node being moved, @Z2 + +
= (cid:11) (cid:12) WZ(cid:22)
@~~ 1@V~2 @x2 (cid:0)" T (cid:0) #
A I: (19) X(cid:0) + (cid:1) +
@x2 ’ 3@x2 @(cid:11) @(cid:12)
(cid:1)(cid:24)W + (cid:1)(cid:17)W : (28)
@x2 @x2
Similarly,thevariationoftheJacobianofthetrans- T (cid:18) (cid:19)
X
f@oWrmZ(cid:22)atio1n@Zsc2ales like one-third of the nodal variation, Asimilarexpressioncanbefgormedforthegyderivative.
@x2 3@x2, andisneglected asa sub-principleterm Neglecting the subsonic blending on M(cid:11) and M(cid:12)
relati(cid:24)ve to the change in solution value, allows (cid:11)+ and (cid:12)+ to be expressed as
@@Wx2j ’WZ(cid:22)@@Zx2j: (20) (cid:11)+ = (cid:11)0;; V~~(cid:11)(cid:11) >00 ; (cid:12)+ = (cid:12)0;; V~~(cid:12)(cid:12) >00 ;
(cid:26) V (cid:20) (cid:26) V (cid:20)
Since the solution is locally assumed to vary only at (29)
@Z1 @Z3
the node being moved, i.e., @x2 = @x2 =0, where ~(cid:11) = V~ n^1 and ~(cid:12) = V~ n^3. The derivatives
V (cid:1) V (cid:1)
3 follow directly,
@Wj @Z2
j=1 @x2 ’WZ(cid:22)@x2: (21) @(cid:11)+ = @@x(cid:11)2; V~(cid:11) >0 ; (30)
X @x2 0; ~(cid:11) 0
@Z2 (cid:26) V (cid:20)
Theremainingtermtoevaluateis @x2. Thesteady-
state distribution of the residual (cid:13)uctuation to the @(cid:12)+ @@x(cid:12)2; ~(cid:12) >0
node can be written = V : (31)
@x2 0; ~(cid:12) 0
(cid:26) V (cid:20)
[(I+M(cid:11))(cid:11) (I +M(cid:12))(cid:12)]WZ(cid:22)Z2 The derivatives of (cid:11) and (cid:12) are
(cid:0)
T
X = [(I +M(cid:11))(cid:11)WZ(cid:22)Z1(cid:0)(I+M(cid:12))(cid:12)WZ(cid:22)Z3]; @(cid:11) = A~y + 1‘1n^1 @A~~ ; (32)
T @x2 2 2 (cid:1)@x2
X
(22) @(cid:11) ~x 1 @~~
= A + ‘1n^1 A; (33)
or @y2 (cid:0) 2 2 (cid:1)@y2
+ + + + @(cid:12) ~y 1 @~~
(cid:11) (cid:12) WZ(cid:22)Z2 = (cid:11) WZ(cid:22)Z1 (cid:12) WZ(cid:22)Z3 ; = A ‘3n^3 A; (34)
T (cid:0) T (cid:0) @x2 2 (cid:0) 2 (cid:1)@x2
X(cid:0) (cid:1) X(cid:0) (2(cid:1)3) @(cid:12) ~x 1 @~~
= A ‘3n^3 A: (35)
@y2 (cid:0) 2 (cid:0) 2 (cid:1)@y2
where
@Z2 @V~2
(cid:11)+ = 1(I +M(cid:11))(cid:11); (24) Having determined @x2, @x2 follows from
2
d(p(cid:26)u) ud(p(cid:26))
(cid:12)+ = 21(I +M(cid:12))(cid:12): (25) du = p(cid:0)(cid:26) ; (36)
The upwinding matrices M(cid:11) and M(cid:12) expressed in and
auxiliary variables reduce to the identity matrix for d(p(cid:26)v) vd(p(cid:26))
supersonic (cid:13)ow, see Ref. 13. Di(cid:11)erentiating Eqn. 23 dv = (cid:0) : (37)
p(cid:26)
while freezing the Jacobians leads to
This completes the derivation of all the terms re-
+ +
@(cid:11) @(cid:12) + + @Z2 quired to evaluate the derivatives of the two-dimen-
WZ(cid:22)Z2+ (cid:11) (cid:12) WZ(cid:22)
@x2 (cid:0) @x2 (cid:0) @x2 sionalobjectivefunction,Eqns.13and14. Reference6
T (cid:20)(cid:18) (cid:19) (cid:21)
X + +(cid:0) (cid:1) presents the analogous derivation in axisymmetric co-
@(cid:11) @(cid:12)
= WZ(cid:22)Z1 WZ(cid:22)Z3 ; (26) ordinates.
@x2 (cid:0) @x2
XT (cid:18) (cid:19) 4of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
4 Ref.6. Thesamereferencealsoshowsthegridconver-
gence of the present solver for this case using uniform
3.5 Free stream mesh re(cid:12)nement, converging in agreement with the
Y, in.
3 benchmark solution.
2.5 Adaption Component Evaluation
2 The individual adaption operations, deletion, swap-
ping, moving,andinsertion,are(cid:12)rst testedseparately
1.5
to gage their relative merits. Each evaluation begins
1 from a converged solution on triangulations of coars-
0.5 ened versions of the structured Laura grid. After
adaptionthesolutionis re-converged. Theadaptation
0 Capsule Sting
will be considered successful if the resultant surface
0 1 2 3 4 5
heat transfer rates are closer to the benchmark data.
X, in. For point deletion the starting mesh is a triangu-
lation of 125 257 structured grid, containing 32,125
Fig.2 ComputationaldomainforMarsPath(cid:12)nder (cid:2)
nodes. The objective for successful point deletion is
capsule.
to remove points without altering the solution. For
Results the (cid:12)rst pass, the adaption is set to remove 10 per-
Case De(cid:12)nition cent (3284) of the nodes, producing minimal change
to the solution. The nodes are predominantly re-
Evaluation of the adaption strategy for hypersonic
movedfromthefreestream. Thenafurther10percent
(cid:13)ow is simpli(cid:12)ed by choosing a two-dimensional lami-
(3157) of the nodes are removed, and the heating
narperfectgascase. Thetestconditionsandgeometry
rates remain essentially unchanged, Figure 3, but the
correspond to the Mach-10 wind tunnel tests of Hol-
15 (cid:13)uctuation splitting scheme has trouble maintaining
lis. The sting-mounted capsule geometry, depicted
the bow shock capture due to the further loss of free
in Figure 2, is modeled as two-dimensional instead
b stream points, as shown in Figure 4. The removal of
of axisymmetric for the present study. The orig-
more points causes rapid solution deteriorationas the
inal axisymmetric wind tunnel model of the Mars
schemefailstoproperly capturethe bowshockonthe
Path(cid:12)ndercapsuleconsistsofaspherically-blunted70-
coarsened meshes.
degree sphere-cone. The body radius is 1 in. and the
nose radius is 21 in. The shoulder radius is 210 in. and For evaluating edge swapping the starting mesh is
the aft body angle is 40 degrees. The base radius is formed from a triangulated 125 129 structured grid
35 in. The sting is 4 in. long with a radius of 1332 in. with 16,125 nodes and 47,868 e(cid:2)dges. Edge swapping
The free stream conditions for the NASA Lang- with the(cid:13)uctuation minimizationstrategyisnot ben-
ley 31-Inch Mach 10 Air Tunnel, corresponding to e(cid:12)cialforthis case. Only about 1percent ofthe edges
a nominal Reynolds number per foot of 0:5 106, are swapped, primarily at the shock, but the result-
are: P = 69 Pa, T = 53 K, (cid:26) = 0:0045 kg/m(cid:2)3, and ing solution does not converge due to ringing of the
u=1416 m/s. The wall temperature is taken to be bow shock near the stagnation point. Figure 5 shows
a uniform 300 K. Historical experience with this tun- pressure contours in the stagnation point region over-
nel at these conditions indicates that laminar perfect laid upon the swapped mesh. Notice the irregular
gascalculationsare adequate for comparisonwith the contours downstream of the shock. The diÆculty in
experimental data. this solution is that the bow shock is making discrete
Themetricforaccuracyischosentobesurfaceheat jumps across grid lines, induced by swapped edges at
transfer rates. The benchmark solution is obtained the shock. The unadapted bow shock more closely
usingthecell-centered(cid:13)uxdi(cid:11)erencesplittingsecond- follows the grid lines at the y = 0 symmetry plane.
order accurateLaura code of Gno(cid:11)o.16 A structured Recall that this mesh was originally well-aligned with
grid with 125 stream-wise nodes and 513 nodes from thebowshockfromthestructured-gridsolver,andthe
thesurfacetotheouterboundary,foratotalof64,125 e(cid:11)ect of edge swapping on a more random unstruc-
nodes, is used to generate the benchmark solution. tured mesh may be di(cid:11)erent. Also, on a (cid:12)ner mesh
The grid is adapted foralignment with the bow shock the smaller gridspacing at the shock might lessen the
and clustered to the boundary layer. The benchmark detrimentale(cid:11)ectsofdiscretejumpsintheshockloca-
solution is shown to be grid converged on this mesh, tion. Theedgeswappingassessmentistriedagainwith
aside from small changes at the stagnation point, in moreconservativethresholds, this timeonlyswapping
halfapercentoftheedges,butstillproducesthesame
b
Reference 6 reports the axisymmetricresults for this case,
disappointing results.
with similarobserved behavior and conclusions as forthe two-
dimensional results, but with the added complexity of special Nodaldisplacementsareevaluatedusingthe16,125-
treatmentfortheaxissingularity. node mesh. The adaption moves 1827 nodes a total
5of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
0.015
0.0125
Y, in.
q
(cid:26)1V130.01 0.5
0.0075
0.005
0.0025 0
-0.5 0
0
0 0.5 1 1.5
X, in.
s, in.
Fig.5 Shockbulgingatstagnationpointafteredge
a)Heatshield.
swaps, pressure contours and grid.
0.001
distanceof 0.7488in., foran averageof 0.0004in. and
RMS of 0.0015 in. This level of movement represents
0.0008
a perturbation to the mesh, with the average move-
(cid:26)1qV013.0006 ment being an order of magnitude smaller than the
edgelengthsattheshock. Themeshmovementoccurs
at the shock, near the shoulder, and in the fore body
0.0004
boundarylayer. Surfaceheatingratesforthiscaseare
shown in Figure 6. The solution is worsened on the
0.0002
fore body while not much change is seen in the wake.
Investigatingthesolutionrevealsthatawindsidevor-
0
tex pattern has emerged, Figure7, that is causingthe
unexpectedheattransferratesontheheatshield. The
-0.0002
1 2 3 4 5 solution is run a further 800,000 iterations with no
changetothiswindsidevortexpattern. Notethatthe
s, in.
solver includes eigenvalue limiting, which is known to
b)Aftbodyandsting. suppressthecarbunclephenomenonforthestructured
grid Roe schemes. The wind side recirculation is not
a physically accurate prediction nor does it appear to
Fig. 3 Surface heating after coarsening,
match the classic carbuncle pathology, but is rather
solid=unadapted, dashed=10% removed, dot-
ted=20% removed. an artifact of discrete mesh point jumps in the bow
shock location, seen in Figure 7 as choppiness in the
bow shock contour lines. For this case the adaption
perturbs a mesh that was initially well aligned to the
4 bow shock in such a way as to largely destroy that
alignment, to the detriment of the solution.
Y, in.
3 Point insertion is tested on the unadapted 16,125-
node mesh. The adaptation adds 3299nodes, approx-
imately a 20 percent increase, with heating results
2
showninFigure8. Manyofthesenewpointsareadded
to the wind side boundary layer. On the heat shield
1
theheatinglevelsrisemarkedlytowardthebenchmark
result, but exhibit a high-frequency oscillation of sig-
0
ni(cid:12)cant amplitude. There is little change in the wake
0 1 2 3 4 5 6
heating.
X, in. Full Adaption Cycle Evaluation
For the full adaption the starting solution is taken
Fig. 4 Loss of shock capture due to coarsening,
onthetriangulated125 129meshwith16,125nodes.
pressure contours. (cid:2)
The intention is to look for an improvement in the
6of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
0.015 0.015
0.0125 0.0125
q q
(cid:26)1V130.01 (cid:26)1V130.01
0.0075 0.0075
0.005 0.005
0.0025 0.0025
0 0
0 0.5 1 1.5 0 0.5 1 1.5
s, in. s, in.
a)Heatshield. a)Heatshield.
0.001 0.001
0.0008 0.0008
q q
(cid:26)1V013.0006 (cid:26)1V013.0006
0.0004 0.0004
0.0002 0.0002
0 0
-0.0002 -0.0002
1 2 3 4 5 1 2 3 4 5
s, in. s, in.
b)Aftbodyandsting. b)Aftbodyandsting.
Fig. 6 Surface heating after nodal displacements, Fig. 8 Surface heating after point insertion,
solid=unadapted, dashed=10% moved, dash-dot- solid=unadapted, dashed=20%inserted, dash-dot-
dot=benchmark. dot=benchmark.
coarsemeshsolutionwithoutincreasingthenumberof
1
nodes. Although not all of the component adaptation
steps were seen to be bene(cid:12)cial, a complete cycle will
Y, in. include each step to allow for synergism between the
components. The strategy for an adaption cycle is to
delete 10 percent of the nodes, swap 1 percent of the
0.5
edges, move 5 percent of the nodes, and then insert
back in 10 percent of the nodes. The solution is re-
convergedbetween each step of the adaption cycle.
The (cid:13)uctuation minimization adaption successfully
removes1601nodes, butrunsintotroubleagainwhile
0 swapping. Windsidevorticesarespawnedinthestag-
-0.5 0 0.5
nationregionbyanoscillatingbowshock. Inane(cid:11)ort
to damp the solution the CFL number is reduced by
X, in. anorderofmagnitude,withoutproducinganimprove-
ment in the solution or eliminating the vortices. The
Fig. 7 Wind Side vortices produced by nodal dis-
number of edges to be swapped is then reduced to
placements, streamlines and pressure contours. 1
2 percent, but the same vortices and oscillating bow
7of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
gas laminar physics as a prototype for the purpose of
0.015
assessing the bene(cid:12)ts of the scheme prior to pursuing
in a three-dimensional reacting gas implementation.
0.0125
The adaption criteria is based upon the minimization
q
(cid:26)1V130.01 of residual (cid:13)uctuations de(cid:12)ned by a multidimensional
upwind distribution scheme for the (cid:13)ow solver. Par-
ticular emphasis is placed on deriving the derivatives
0.0075
of the nodal displacement objective function for the
upwind (cid:13)uctuation splitting distribution algorithm.
0.005
The hypersonic demonstration case is a two-dimen-
sionalsting-mountedcapsuleatlaminarMach10wind
0.0025
tunnel conditions, modeled with perfect gas air. Sur-
face heat transfer is the quantity of interest used to
0
0 0.5 1 1.5 measure the success of the solution adaptive process.
Four elementary operations are performed during
s, in.
the adaption: point deletion, edge swapping, nodal
a)Heatshield. displacement, and point insertion. The limits of ef-
fectiveness for each of the four operations is probed
0.001
independently to guide their relative weighting in the
full adaption strategy. Coarsening of the mesh works
0.0008
wellremovinguptoabout20percentofthenodes,but
(cid:26)1qV013.0006 leadsto aloss of the bow shockcapture and a rapidly
deteriorating solution for more aggressive coarsening.
The other three adaption techniques generally have
0.0004
minimal or negative impacts on the solution, with
the exception of point insertion where solution im-
0.0002
provement occurs but at the cost of high frequency
oscillations. Overall a lack of robustness is demon-
0
strated,usuallycausedbydistortionstothebowshock
nearthe stagnationpoint. Common failuremodes are
-0.0002 an oscillating shock that sheds wind side vortices, or
1 2 3 4 5
a steady shock kink that produces unexpected pres-
s, in.
sure contours and surface heating spikes. For the full
b)Aftbodyandsting. adaptioncycletheforebody heatingimprovesbutbe-
comesoscillatory,whiletheaftbodyandstingheating
predictions improve slightly.
Fig. 9 Surface heating after full adaption cy-
One likely cause of the lack of robustness exhibited
cle, solid=unadapted, dashed=adapted, dash-dot-
dot=benchmark. by the adaption scheme at the stagnation point bow
shock is the con(cid:13)ict between characteristic directions
shock appear. To proceed with this case the edge between the supersonic upstream (cid:13)ow, which is per-
swapping is omitted entirely and the adaption cycle fectly upwinded by the (cid:13)ow solver, and the subsonic
continued with the nodal displacement step, where shocklayer,wheretheacousticcharacteristicsapprox-
706 nodes are moved a total distance of 0.7 in. Fi- imated asdiscretewaves. A second contributorto the
nally, 1638 nodes are added, yielding the results of lack of robustness is the local scope of the adaption
Figure 9. The fore body heating is generally im- strategy. Theschemeisintentionallylocaltobesyner-
provedtowardthebenchmarksolution,althoughthere gisticwiththe(cid:13)owsolverwhichhascompactsupport,
is a high-frequency oscillation in the data starting at sharing the same software data structures and loop
s = 0:6. The aft body heating is improved to match controls. The local scope also is desirable for paral-
thebenchmarkbetweens=1:2{1.8,althoughtheorig- lel processing. However, a larger stencil, an implicit
inalsolutionwasnotfaro(cid:11)fromthebenchmarkinthis scheme, or an adjoint equation criteria may be neces-
region. The sting heating is only slightly changed aft sary to provide suÆcient smoothness at shocks.
of s=3.
Although heat transfer is a variable of primary in-
terest for aerothermodynamics and may be a (cid:12)rst
Concluding Remarks consideration for driving the adaption criteria, the
Ananisotropicunstructuredmeshadaptionstrategy present results emphasize that for blunt-body hyper-
for hypersonic viscous (cid:13)uid (cid:13)ows is presented. The sonic aerothermodynamics both the bow shock loca-
strategy is implemented for two-dimensional perfect tionandsmoothnessarecriticaltothecomputationof
8of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
13
the forebody heat transferrates. Due to the strength Wood, W. A. and Kleb, W. L., \2-D/Axisymmetric For-
mulation of Multi-dimensional Upwind Scheme," AIAA Paper
of the hypersonic bow shock, discrete jumps in shock
2001-2630,June2001.
position from one node to another produce signi(cid:12)cant 14
Wood, W. A., Multi-dimensional Upwind Fluctuation
pressuredisturbancewavesandentropyvariationsthat Splitting Scheme with Mesh Adaption for Hypersonic Viscous
stronglya(cid:11)ectthe(cid:13)ow(cid:12)eldinthesubsonicbubbleand Flow,Ph.D.thesis,VirginiaTech,2001.
15
overwhelmmore subtle adaptionwithin the boundary Hollis, B. R. and Perkins, J. N., \Transition E(cid:11)ects on
Heating in the Wake of a Blunt Body," Journal of Spacecraft
layer.
and Rockets,Vol.36,No.5,Sept. 1999,pp.668{674.
To summarize, the present adaption techniques do 16
Gno(cid:11)o, P.A.,Gupta, R.N.,andShinn,J. L.,\Conserva-
not signi(cid:12)cantly improve the solutions, the adaption tionEquationsandPhysicalModelsforHypersonicAirFlowsin
e(cid:11)ectiveness is not consistent, and the solutions dis- ThermalandChemicalNonequilibrium,"NASATP2867,Feb.
1989.
play a lack of robustness. Point deletion works the
best, but is risky if applied too aggressively. Uniform
re(cid:12)nement is both simpler and more reliable than the
solution adaptive strategy embodied by the present
method.
References
1
Thomas,J.L.,Alexandrov,N.,Alter,S.J.,Atkins,H.L.,
Bey,K.S.,Bibb,K.L.,Biedron,R.T.,Carpenter,M.H.,Cheat-
wood, F. M., Drummond, P. J., Gno(cid:11)o, P. A., Jones, W. T.,
Kleb,W.L.,Lee-Rausch, E.M.,Merski,N.R.,Mineck,R. E.,
Nielsen, E. J., Park, M. A., Pirzadeh, S. Z., Roberts, T. W.,
Samareh, J. A., Swanson, R. C., Vatsa, V. N., Weilmuenster,
K.J.,White,J.A.,Wood,W.A.,andYip,L.P.,\Opportuni-
tiesforBreakthroughsinLarge-ScaleComputationalSimulation
andDesign,"NASA/TM2002-211747,June2002.
2
Venditti,D.A.andDarmofal,D.L.,\GridAdaptationfor
Functional Outputs: Application to Two-Dimensional Inviscid
Flows,"JournalofComputationalPhysics,Vol.176,No.1,Feb.
2002,pp.40{69.
3
Habashi, W.G., Fortin, M.,Dompierre,J., Vallet, M.-G.,
and Bourgault, Y., \Anisotropic Mesh Adaption: A Step To-
wards a Mesh-Independent and User-Independent CFD," Bar-
riersandChallengesinComputational FluidDynamics,edited
byV.Venkatakrishnanetal.,KluwerAcademicPublisher,1998,
pp.99{117.
4
Habashi,W.G.,Dompierre,J.,Bourgault,Y.,Fortin,M.,
andVallet,M.-G.,\Certi(cid:12)ableComputationalFluidDynamics
Through Mesh Optimization," AIAA Journal, Vol. 36, No. 5,
May1998,pp.703{711.
5
Wood, W. A.and Kleb, W. L., \Di(cid:11)usion Characteristics
of Finite Volume and Fluctuation Splitting Schemes," Journal
ofComputational Physics,Vol.153,No.2,Aug.1999,pp.353{
377.
6
Wood, W. A., \Multi-dimensional Upwind Fluctuation
Splitting Scheme with Mesh Adaption for Hypersonic Viscous
Flow,"NASA/TP2002-211640,April2002.
7
Caraeni,D.,Caraeni,M.,andFuchs,L.,\AParallelMul-
tidimensional Upwind Algorithm for LES," AIAA 2001-2547,
June2001.
8
Roe, P., \Compounded of Many Simples," Barriers
and Challenges in Computational Fluid Dynamics, edited by
V. Venkatakrishnan et al., Kluwer Academic Publishers, 1998,
pp.241{258.
9
Roe,P.andNishikawa,H.,\AdaptiveGridGenerationby
MinimisingResiduals," ICFD Conference on Numerical Meth-
ods in Fluid Dynamics, Oxford, U.K., 2001, submitted to the
InternationalJournalforNumericalMethodsinFluids.
10
Wood, W. A. and Kleb, W. L., \On Multi-dimensional
UstructuredMeshAdaption,"AIAAPaper99{3254,June1999.
11
Yamaleev, N. K. and Carpenter, M. H., \On Accuracy
of Adaptive Grid Methods forCaptured Shocks," submitted to
Journal of Computational Physics,2001.
12
Sidilkover, D., \Multidimensional Upwinding and Multi-
grid,"AIAAPaper95{1759,June1995.
9of9
AmericanInstitute ofAeronauticsandAstronauticsPaper2003{3824
