Adjoint-Based, Three-Dimensional Error Prediction and Grid Adaptation

12 Pages·2002·1.1 MB·English

by  Park, Michael A

AIAA 2002-3286
Adjoint-Based, Three-Dimensional
Error Prediction and Grid Adaptation
Michael A. Park
NASA Langley Research Center, Hampton, VA 23681
32nd Fluid Dynamics Conference
24-27 June 2002
St. Louis, MO
For permission to copy orrepublish, contact the American Institute ofAeronautics andAstronautics
1801Alexander Bell Drive, Suite 500, Reston, VA20191-4344
Adjoint-Based, Three-Dimensional Error
Prediction and Grid Adaptation
Michael A. Park *
NASA Langley Research Center, Hampton, VA 23681
Abstract (3D) configurations compared with structured grid
Engineering computational fluid dynamics (CFD) techniques. The use of unstructured grid technology
analysis and design applications focus on output func- for CFD simulations allows more freedom in adapt-
tions (e.g., lift, drag). Errors in these output functions ing the discretization of the meshes to improve the
are generally unknown and conservatively accurate so- fidelity of the simulation. Many previous efforts have
lutions may be computed. Computable error estimates attempted to tailor the discretizations of unstructured
can offer the possibility to minimize computational meshes to increase solution accuracy while reducing
computational cost. 18
work for a prescribed error tolerance. Such an esti-
mate can be computed by solving the flow equations Most of these adaptive methods focus on modifying
and the linear adjoint problem for the functional of discretizations to reduce local equation errors. These
interest. The computational mesh can be modified local errors are not guaranteed to directly impact er-
to minimize the uncertainty of a computed error es- ror in global output functions. These methods, often
timate. This robust mesh-adaptation procedure au- referred to as feature-based adaptation, focus on re-
tomatically terminates when the simulation is within solving discontinuities or strong gradients in the flow
a user specified error tolerance. This procedure for field. Flow features (e.g., shocks) can be in the in-
estimating and adapting to error in a functional is correct location due to errors elsewhere in the flow
demonstrated for three-dimensional Euler problems. field. Also, resolving the flow in a location may have a
An adaptive mesh procedure that links to a Computer minimal effect upon the output function (e.g., a down-
Aided Design (CAD) surface representation is demon- stream shock).
strated for wing, wing-body, and extruded high lift Ifthe flow equations are linearized about the flow so-
airfoil configurations. The error estimation and adap- lution, the solution of a linear dual problem can yield a
tation procedure yielded corrected functions that are direct measure of the impact of local primal (flow equa-
as accurate as functions calculated on uniformly re- tion) residual on a selected functional output. The
fined grids with ten times as many grid points. combination of the primal and dual problems can also
be applied to yield a correction to a specified func-
Introduction tional on a given mesh. If a specified error tolerance
in an output function is required, the cost of comput-
Engineering problems commonly require computa-
ing a CFD solution can be minimized by adapting the
tional fluid dynamics (CFD) solutions with functional
discretization of the problem to directly minimize un-
outputs of specified accuracy. The computational re-
certainties in the corrected output function of interest.
sources available for these solutions are often limited
Also, the entire adaptive simulation can be terminated
and errors in solutions and outputs are often unknown.
when the predicted error is equal to a specified toler-
CFD solutions may be computed with an unnecessarily
ance, preventing the waste of computational work on
large number of grid points (and associated high cost)
an overly large mesh.
to ensure that the outputs are computed to within a
There are many examples of these techniques in the
required accuracy. A method to estimate the error
finite element communities. °'10 Pierce and Giles 11ap-
present in a computed functional offers the possibility
plied these methods to finite-volume discretizations.
to avoid the use of overly refined grids to guarantee
Venditti and Darmofa112'13 have demonstrated these
accuracy.
methods for compressible two-dimensional (2D) invis-
Unstructured grid technology promises easier initial
cid and viscous flow solutions. The present study is
grid generation for novel complex three-dimensional
essentially a 3D extension of the methods of Venditti
and Darmofa113 that adapts the mesh to reduce uncer-
*Research Scientist, Computational Modeling and Simula-
tion Branch, m. a. park@larc, nasa. gov tainty in an error correction. Miiller and Giles 14have
Copyright @ 2002 by the American Institute of Aeronautics and also presented results for a similar approach utilizing
Astronautics, Inc. No copyright is asserted in the United States
under Title 17, U.S. Code. The U.S. Government has a royalty- a different adaptation parameter that directly targets
free license to exercise all rights under the copyright claimed herein the error correction.
for Governmental Purposes. All other rights are reserved by the
copyright owner. The use of adjoint variables (solution of the dual
1 OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
problem) is an efficient method for computing deriva- Adjoint Equations
tives of a functional of interest for gradient-based After the flow solution is known, the discrete ad-
design methods. Some examples of discrete adjoint de- joint equations 15'is, 19are solved to complete the dual
sign methods are given in Nielsen and Anderson. 15'16 problem. The first step isto linearize the flow equation
The discrete dual problem for adjoint variables can be residual R and output function f with respect to the
expensive. However, the adjoint solution is already flow solution Q. After this linearization, an adjoint
available for use during the aerodynamic design pro- variable A is solved for each of the flow equations.
cess, so it could be employed for simultaneous design, An abbreviated derivation, adapted from Taylor et
error prediction, and grid adaptation. al., 22is below. The chain rule for the linearized output
The combination of adjoint-based grid adaptation function is
and design techniques can yield an attractive tool (0 ' TOR
for the aerodynamic design of new configurations.
Adjoint-based error prediction and adaptation can
The adjoint variable A is defined as the effect of the
yield smaller meshes than traditional feature-based
flow residual on the output function:
schemes with computable error estimates on out-
put functions. Design processes require analysis and
Of A (4)
derivative evaluation tools that operate with mini- OR
real human interaction. Robust, automatic adaptation
techniques enable the increased use of nonlinear flow A set of linear equations is solved to find A:
calculations in larger multidisciplinary design frame-
works. These new techniques will enable efficient
oq/ (5)
analysis for existing configurations and expanded ex-
ploration of design spaces for new configurations.
After the flow solution isknown, this set of linear equa-
tions is solved with GMRES. 23See Refs. 15,16, and 19
Flow Equations for details. A implicit point-iterative time-marching
The FUN3D lr 19t (Fully Unstructured Navier- method is employed to compute the adjoint solution
Stokes Three-Dimensional) suite of codes is employed for the high lift configuration. 2426
in this study. The compressible flow solver employs
an unstructured finite-volume tetrahedral method for Error Correction
conserved variables, Q, i.e., The error prediction and correction scheme is taken
from Ref. 13. With a solution on a mesh of reasonable
Q [p pu pv pw E]T (1) size Q0, it is desirable to predict the value of an out-
put function evaluated with a solution on a much finer
where p is density, u, v, and w are velocity, and E is mesh f(Q*). This prediction can be computed with-
total energy per unit volume. The incompressible flow out the solution on this finer mesh when the adjoint
solver employs the following state variables: solution on this reasonable mesh A° is utilized. The
full derivation of the error correction term is available
in Refs. 11, 13, 14. An abbreviated derivation is pre-
sented by expanding the Taylor series about f(Q0),
where p is pressure. The node-based variables Q are i.e.,
computed by driving the flow equation residual R to
steady-state with an implicit point-iterative method.
The code is able to solve incompressible, Euler, and /(q.) /(qO)+ _ 0 "'"
Reynolds-Averaged Navier-Stokes (RANS) flow equa- (6)
tions loosely coupled to the Spalart-Allmaras 2° one- Employing eq. (4) and assuming that the residual
equation turbulence model. The present study em- on the much finer mesh is zero yields an improved es-
ploys only the Euler and incompressible equations. timate for the functional of interest f_t:
The solution of Q allows the calculation of integral
quantities f (e.g., lift, drag). To speed execution, the
problem domain is decomposed and the flow and the OORf o A° (7)
adjoint problems are solved with a parallel execution
scheme utilizing the Message Passing Interface (MPI) R(Q*) 0 (8)
standard. The FUN3D suite of codes is being extended f(q*) v f_t f(q0)_ (),0)rR(q0) (9)
to the HEFSS (High Energy Flow Solver Synthesis) 2s
modular framework of FORTRAN 90 shared libraries. To improve the prediction of the functional f_t, Q0
and A° can be interpolated to an embedded mesh. In-
thttp://fun3d, larc.nasa.gov terpolation is performed in two ways for this study:
2 OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
a linear interpolation (QL,,_L) and a higher order in- The intensity I is therefore the average of a primal
terpolation (QH,AH). Details of this interpolation and residual weighted with a dual solution interpolation
the construction of this embedded mesh are in the "In- error and a dual problem residual weighted with a
terpolation Techniques" and "Embedded Mesh" sec- primal solution interpolation error. This form of the
tions. Substituting these interpolated quantities into adaptation intensity tends to focus on the nonlinear
eq. (9) yields the linear and higher order functional contributions to the function error, which increases ro-
estimates fL and f_t: bustness of the adaptation method.
f c /(QC)_ (AL)TR(QC) (10)
Error Correction and Adaptation
fHt f(QH) _ (AH)TR(QH) (11) Process
The error correction and adaptation process begins
Adaptation Parameter
with an initial tetrahedral mesh, which can come from
The adaptation parameter, also from Ref. 13, is in-
any mesh generation system. The state variables are
tended to specify a grid spacing modification to reduce
computed as the nonlinear solution to the flow equa-
the uncertainty in the computed error prediction. The
tions on the initial mesh. The adjoint variables are
grid is not modified to directly reduce the computed
then computed with the linearized flow equations at
error prediction (as in Ref. 14) because an estimate for
the flow solution. These flow and adjoint solution pro-
the functional including this error term can be com-
cedures employ a parallel execution scheme. Then the
puted with eq. (9). Instead, targeting the uncertainty
global problem domain is reconstructed to facilitate
in this computed quantity is more effective and im-
the creation of a finer, embedded grid with interpo-
proves the robustness of the adaptive process. The
lated primal and dual solutions.
error correction (eq. (9)) including the uncertainty in
the dual solution is Embedded Mesh
To compute the error prediction and the adapta-
f(Q0) _ f(Q,) v (A0)TR(Q0) + (A* -- A°)TR(Q °)
(12) tion parameter a globally embedded or h-refined mesh
is created. To construct the embedded mesh, a new
The uncertainty in the computed error correction is
node (open circle) is inserted at the midpoint of each
f_t - f(Q*) v (A* - A°)TR(Q °) (13) existing edge that connects existing nodes (closed cir-
cles); see Fig. l(a) and l(b). Each existing tetrahedron
The relation of the primal and dual problems 11,13 is subdivided to reconnect these new nodes with eight
yields another expression for the error correction un- interior tetrahedra. (Each of the existing boundary
certainty faces is also divided into four triangles.)
(a* - a°)rR(q °) R_(a0)(q, _ q0)r (14)
Where Rx(A) is the residual of the dual problem:
(0R T_ (15)
\OQ] \OQ]
A computable term is found by using the interpolation
error of A to replace (A* - A°) and the interpolation
error of Q to replace (Q* - Q0). The higher order
a) Original tetrahedron. b) Embedded tetrahe-
interpolate for Q0 and A° is employed to improve pre- dron.
diction in place of the linear interpolate in Ref. 13. Fig. 1 Tetrahedron embedding process.
The interpolation error is expressed as the difference
in the high-order and linear interpolated values: The four new tetrahedra constructed at the corners
of the existing tetrahedron have the same shape as
(a* - _,o) _ (_,H _ _,L) (16)
the original tetrahedron but are smaller in size. The
(q. _ q0) _ (QH_ QL) (17) construction of the four corner tetrahedra leaves an
interior volume with eight faces, which is subdivided
The average of the absolute values of the two uncer- into four tetrahedra. The four interior tetrahedra have
tainty terms in eq. (14) yields the adaptation intensity three unique configurations. The configuration with
I, which is computed for each equation on each em- the lowest maximum dihedral angle is selected.
bedded node:
The new nodes are placed at the midpoints of edges
during the mesh embedding process. The embedded
I I(AH -- AL)TR(QH)I + IRx(AH)(QH - QL)TI
nodes on the boundaries of the mesh may no longer
(18) remain on the original surface definition of the model.
3OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
When the grid is adapted to improve the discretiza- Error Correction and Adaptation Parameter
tion, the surface fidelity of the mesh is maintained with Once the new embedded grid is constructed with
boundary node projection. QL, AL, QH, and AH, it is partitioned to allow paral-
Interpolation Techniques lel calculation of the functional and flow and adjoint
equation residuals. The flow and the adjoint equations
The primal and dual solution variables y, which are
are not iterated or solved on this embedded grid; the
all elements of Q and A, are now interpolated to this
embedded mesh. The value of the solution at each flow state and adjoint variables are interpolated from
the original mesh. Therefore, the only computational
of the existing nodes is directly copied into the corre-
costs on this larger embedded grid are function evalu-
sponding nodes of the embedded mesh. Each of the
ations, flow and adjoint residual evaluations, and dot
solution variables is interpolated in two ways to form
products of vectors. The linear and higher order er-
the linear and the higher order reconstruction for the
ror correction term, eq. (10) and (11), is computed at
new nodes. The higher order reconstruction of the so-
each node on the embedded mesh and summed over
lution for the new nodes requires the computation of
the entire mesh for all flow equations.
least-squares gradients 17 at the existing nodes using
The adaptation intensity, eq. (18), is also computed
the existing mesh and solution. To simplify the 3D
at each embedded mesh node. At each node that is
interpolation, the interpolation is performed by inde-
also present in the original mesh, the computed inten-
pendently examining each existing edge in the original
sity is zero due to the chosen interpolation technique.
mesh; thus the interpolation problem becomes one-
The values of the lower and higher order interpolation
dimensional along each existing edge. Each existing
schemes are the same at these existing nodes; thus
edge has an existing node at each edge endpoint and
a new node at the edge midpoint (see Fig. 2). (QH _ QL) and (,_H _ ,_L) are exactly zero.
To specify the grid adaptation on the original mesh,
_____-0
the adaptation intensities must be reduced from the
_x2, 5'2,_)2 embedded mesh to the original mesh (I0). The new
nodes on the embedded mesh all lie on existing edges
Y3/2
of the original mesh (see Fig. l(b)). Therefore, to con-
xz, Yz, Yz
struct I0 the original mesh is examined one edge at
a time (see Fig. 2). One half of the intensity com-
Fig. 2 New node and existing edge. puted at each new node (which is at the midpoint of
these original edges) is added to each existing node at
An edge has two 3D endpoints: xl and x2. The
the endpoints of these edges. The intensities are also
vector that represents the length and direction of the
summed over the equations at this point, resulting in
edge is Ax x2-xl. The 3D least-squares gradient of
one intensity value for each original node.
the solution Ay can be projected to a total derivative
The adaptation parameter, which has been reduced
along the edge to facilitate interpolation by
to the original mesh, is summed to find the global
_I VYTAX (19) intensity Iy _ I0. The number of nodes in the orig-
inal mesh n and the user-specified error tolerance t are
The new node interpolation Y3/2 can be expressed combined to scale the adaptation intensity; that is
as a combination of the solution values and the deriva-
tives at each endpoint (Yl, _)1, Y2, and _)2). The linear Ig_ (22)
interpolation Y3L/2 is the average of the two end nodes I_ TTI0
To perform the grid adaptation, the mesh is locally
L Yz + Y2 (20)
Y3/2 2 enriched in the location of nodes where the scaled
intensity I_ is greater than a value, e.i., unity. Ref-
The higher order interpolation Y_2 is found with a
erence 13 demonstrates a way to specify a new el-
cubic fit of the endpoint data that is evaluated at the
ement spacing function. The adaptation procedure
midpoint
self-terminates as all elements of I_ become less than
H Yz + Y2 _)1- _)2 (21) unity (i.e. no nodes are flagged for adaptation).
Y3/2 2 _ 8
Equation (21) is equivalent to aleast-squares quadratic
fit of the endpoint data that is evaluated at the edge Adaptation Mechanics
midpoint. The adaptation mechanics utilize three independent
At the completion of the grid embedding and inter- modules. The first module inserts new nodes into the
polation step, the linear and higher order interpolated existing mesh and locally reconnects tetrahedra and
solutions to the primal and dual problems QL, AL, boundary faces to maintain a valid tessellation. The
QH, and AH are available to compute the error cor- second module employs face and edge swapping to im-
rection and adaptation parameter. prove the mesh quality. The final module performs
4OF ii
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
grid smoothing and boundalT/node projection opera- plemented by Brasher and Park. 31 This initial imple-
tions. mentation has been extended and coupled to CAPRI,
which utilizes native CAD point projection routines.
Node Insertion
Nodes on the boundary are smoothed by moving the
The node insertion method is currently one level
nodes in CAD (u, v) parametric space to improve the
of selective h-refinement. To start the refinement, all
shape measure of adjacent tetrahedra.
the edges surrounding nodes on the original mesh that
As the nodes are projected, the neighboring tetra-
have ascaled intensity I_ greater than unity are flagged
hedra are tested for validity. If invalid tetrahedra
for h-refinement. The set of flagged edges is exam-
resulted from the projection, the projection distance
ined tetrahedron by tetrahedron and additional edges
of the boundary nodes is reduced until the neighboring
are flagged in an attempt to maintain grid quality
tetrahedra are valid. Then the nodes in the neighbor-
(i.e., low maximum dihedral angles, few high-degree
hood of the projected node are smoothed to improve
nodes). The final set of flagged edges results in tetra-
a quality measure of the adjacent tetrahedra. The
hedra with one edge, three edges on one face, or all six
boundary points are then moved into the fully pro-
edges flagged. A tetrahedron with all six edges flagged
jected position in a number of iterative cycles.
is illustrated in Fig. l(b). The mesh is then h-refined
It is anticipated that grid smoothing in the neigh-
by inserting new nodes on the midpoint of the flagged
borhood of projected nodes may not adequately re-
edges and reconnecting these nodes into new tetrahe-
gain surface fidelity of anisotropic meshes. A grid-
dra and boundary faces.
movement scheme may be required as in Ref. 16. An-
Face and Edge Swapping other possibility is a 3D version of mesh restructuring
as in Ref. 32.
The current post-adaptation grid-improvement
scheme employs face and edge swapping. 27 The swap- Adaptation Module Interaction
ping algorithm minimizes a shape (cost) function (e.g.,
The current selective h-refinement technique often
aspect ratio or dihedral angle). This study sought to creates high-degree nodes on the border of the adapted
reduce only the cell aspect ratio AR
regions. The smoothing algorithm is currently unable
to improve elements that are adjacent to high-degree
i tetrahedral eireumsphere radius
AR (23) nodes. The edge and face swapping techniques effec-
3 tetrahedral in-sphere radius
tively improve shape measures and reduce the number
Reconnections of tetrahedra with undesirable shape of high-degree nodes, facilitating projection and node
measures are investigated and new local tetrahedra smoothing.
configurations with more desirable shape measures These three adaptation modules where developed
are selected. Edges on boundary faces can also be independently to facilitate a quick initial implementa-
swapped. To simplify and speed up the edge swapping tion and to investigate the strengths and weaknesses
routine, the boundary face information is discarded of each technique. Merging the abilities of these three
and reconstructed at the end of the swapping process. separate modules will allow for more flexible modifica-
Smart-Laplacian smoothing 2s is used on the interior tions of grids (e.g., point insertion, point removal, and
nodes. The actual locations of the boundary nodes is anisotropic elements).4, 6
not modified in this module; that modification is per-
Results
formed by the grid smoothing and projection module.
Adaptation results are shown for a wing, wing-body,
Grid Smoothing and Projection
and high lift configurations. The wing is simulated
The inserted boundary nodes may not be located
with incompressible and transonic flow conditions.
on the surface geometry of the model to be simulated
The wing-body and high lift configurations are sim-
since they were inserted at the midpoints of existing ulated with subsonic flow.
edges. A CAD model is employed to describe the
actual model surface. To regain the surface fidelity Initial Mesh Generator
of the mesh, the newly inserted boundary nodes are The initial meshes for these error prediction and
projected to the model surface with a CAD interface adaptation studies are generated with the FELISA
package CAPRI (Computational Analysis PRogram- mesher 4 connected to CAD geometxT/ by CAPRI
ing Interface).29'3° The projection of these new nodes through the GridEx 33framework. FELISA is a Delau-
to their location on the CAD surface can result in in- nay mesh generator with an advancing-front method
verted, invalid tetrahedral elements. for inserting nodes. The GridEx framework is cur-
A grid node-smoothing algorithm is employed to rently being developed at NASA Langley Research
facilitate boundary projection without generating in- Center to link various grid generation and adaptation
valid elements and to improve the overall quality strategies to geometry through CAPRI. This frame-
(shape measure) of the mesh. This post-adaptation work is also utilized in a batch mode to perform uni-
smoothing is similar to F_'eitag 2s and was initially ira- form grid refinement studies.
5 OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
0.1 r r
c_ Second-Order Slope
0.01
rm
\\
E
o User-Specified Error Tolerance ',_\
o 0.001
',@
'a
DragAdaptation,Tol0.001_ ] _>
WithLinearCorrection _ /
WithCubicCorrection -o
0.0001
1000 10000 100000 le+06
Number of Nodes
Fig. 4 Coefficient of drag for the ONERA M6
a) Initial ONERA M6 mesh. adapted to incompressible drag.
gle of attack of 0 deg. Directly computed drag and
estimates of drag for an ONERA M6 wing as a func-
tion of number of nodes is shown in the Fig. 4 log-log
plot. The adaptation and error correction results are
shown for a drag error tolerance of 0.001. The directly
computed drag on the adapted meshes is represented
by the solid lines with circular symbols. The error-
corrected drag calculated with the linear interpolated
solution f_t is represented by a dashed line and square
symbols. The estimated functional calculated with the
higher order interpolated solution f_t is represented
by a dashed line and diamond symbols. The correct
drag is zero because of non-lifting, subcritical, invis-
cid flow. Therefore, the y-axis denoted "Coefficient of
Drag" is also the error in drag.
The triangle labeled "Second-Order Slope" in Fig. 4
illustrates second-order spatial convergence. The
adaptive-grid method results in drag calculations
b) ONERA M6 mesh adapted to incompressible drag. that converge at a much higher rate than the
asymptotic convergence rate of uniform refinement
Fig. 3 Initial and adapted ONERA M6 meshes.
(second-order). The adaptive procedure correctly self-
Drag Adaptation - Incompressible ONERA M6 terminated when the drag error of the adapted flow
solution reached the user specified error tolerance (dot-
The initial mesh for an ONERA (Office National
dash line).
d'Etudes et de Recherches A6rospatiales) 5/[6 wing
with 5227 nodes is shown in Fig. 3(a). The mesh The final grid (454 thousand nodes) after five cycles
has extremely coarse spacing, especially at the trailing of grid adaptation to incompressible drag is shown in
edge and is intended to resolve the surface curvature of Fig. 3(b). The adaptation process clustered grid points
the leading edge and the wingtip and not any particu- at the leading and trailing edges of the wing. Points
lar flow features. The spacing function for this mesh is are also clustered in the neighborhood of the stagna-
specified manually and is intended to be representative tion stream line.
of an automated curvature or maximum chord-height The tetrahedra shape measure AR (eq. 23) is min-
specification. The CAD geometry is represented with imized by the mesh improvement algorithm. The
the CAPRI native kernel FELISA with a part gener- boundary node smoothing algorithm is intended to op-
ated fl'om a surface IGES definition. timize the shape measures of the tetrahedral elements.
The initial ONERA 5/[6 mesh was used in the grid Therefore, the shape measures of the boundary faces
adaptation process with incompressible flow at an an- depicted in this surface plot may not be optimal.
6 OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
0.035 0.035
] [
Uniform Refinement Uniform Refinement
With Linear Correction [] Drag Adaptation
With Cubic Correction o With Cubic Correction <>
0.03 0.03 Extrapolated
c_\\\ Extra olated
D \
0.025 0.025
rm rm
'b _6
"_ 0.02 "_ 0.02
oo 0.015 0o 0.015
0.01 0.01
0.005 0.005
1000 10000 100000 le+06 1000 10000 100000 le+06
Number of Nodes Number of Nodes
a) Uniformly refined drag convergence. a) One level of h-refinement.
0.018 0.035 ] [
Uniform Refinement
With Linear Correction [] I Drag Adaptation
0.017 With Cubic Correction <> I With Cubic Correction _>
0.03 Extrapolated
0.016
0.025
0.015
C3 _/5///
rm
/ 'b
0.014 E 0.02
._o
0.013
oo o
o 0.015
0.012
0.011 0.01
0.01 I I I I I I I I I
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.005
1000 10000 100000 le+06
(h0/h) 2
Number ofNodes
b) Richardson extrapolation. b) Two levels of h-refinement.
Fig. 5 Coefficient of drag for the uniformly refined Fig. 6 Coefficient of drag for the drag adapted
ONERA M6 at 0.84 Mach. ONERA M6 at 0.84 Mach.
Drag Adaptation - Transonic ONERA M6 finement of the ONERA M6 wing. A representative
A uniform refinement of the ONERA M6 wing mesh element length h was estimated by computing the cube
is computed at 0.84 Mach and an angle of attack of 3 root of the number of nodes. This length was nor-
deg. The drag directly computed by the flow solver realized by the estimated length of the 624 thousand
is shown with linear and higher order interpolated er- node mesh h0. The symbols are drag computed by
ror corrections as a function of the number of nodes the flow solver and error corrected values. A linear
in Fig. 5(a). The extrapolated (grid-converged) drag fit of the data at (ho/h) 2 1.0 and (ho/h) 2 1.7
value for Fig. 5(a) was estimated with Richardson ex- is used to estimate the grid-converged answer for all
trapolation from Fig. 5(b). These meshes have the three schemes. All three schemes indicate a similar
same spacing function as Fig. 3(a) globally modified grid-resolved value. An additional flow solution (1.2
with a scalar to uniformly reduce the element spac- million nodes, (ho/h) 2 0.6) is shown to verify the
ing. These grids were generated with the batch version accuracy of the linear fit; it was not used to construct
of GridEx using the FELISA mesher and CAPRI for the linear fit of the computed drag. An improved in-
CAD geometry access. terpolation scheme for the error correction would yield
Figure 5(b) shows drag and estimates of drag as a superconvergent functional estimate as in Ref. 11.
a function of element size for the uniform grid re- The initial coarse mesh shown in Fig. 3(a) is era-
7OF ii
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
a) Upper wing surface mesh. a) Upper wing surface mesh.
b) Upper wing surface Mach contours. b) Upper wing surface Mach contours.
Fig. 7 Initial ONERA M6 upper wing surface. Fig. 8 Final ONERA M6 upper wing surface,
adapted with one level of h-refinement.
ployed in the error prediction and grid adaptation
procedure with two different adaptation methods. The tion convergence is illustrating the limitations of using
coefficient of drag is plotted as a function of mesh size a single level of selective h-refinement as the adap-
in Fig. 6. The user specified error tolerance in drag tive node-insertion procedure. The use of two levels of
is 0.0019. The unifornfly refined flow solution from h-refinement better approximates a continuous vari-
Fig. 5(a) is shown with the adapted grid flow solu- ation in element size. The two adaptation methods
tion and higher order error prediction of the adapted converged to similar meshes and drag values.
grid in Fig. 6. Figure 6(a) demonstrates a single level The upper wing surface grid and Mach contours
of h-refinement for all nodes with a scaled adaptation of the initial flow field computed on the mesh from
intensity I_ greater than one. Figure 6(b) shows h- Fig. 3(a) is shown in Fig. 7. The shocks in this
refinement for I_ greater than one and a recursive call initial grid are poorly resolved. The mesh (357 thou-
to the adaptation mechanics for I_ greater than 75, sand nodes) and Mach contours of the ONERA M6
yielding two levels of h-refinement at each adaptation adapted to drag with one level of selective h-refinement
cycle. The initial convergence of the function is better is shown in Fig. 8. The mesh (374 thousand nodes)
in Fig. 6(b) than Fig. 6(a). This improvement in func- and Mach contours of the ONERA M6 adapted to
8OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286
a) Upper wing surface mesh. Fig. 10 Initial EET grid.
UniformRefiInement
0.44 LiftAdaptation
WithLinearCorrection <>
WithCubic Correction
Y
X 0.42
0.4
s--
o J
O z
0.38
S
0.36
10000 100000 le+06 le+07
Number ofNodes
Fig. 11 Coefficient of lift for the lift-adapted EET
b) Upper wing surface Mach contours.
at 0.40 Mach.
Fig. 9 Final ONERA M6 upper wing surface,
adapted with two levels of h-refinement. surface details of the fuselage, wing leading edge, blunt
trailing edge, and the wingtip. The geometry is repre-
drag with two levels of selective h-refinement is shown sented with a Parasolid CAD kernel accessed though
in Fig. 8. The final meshes and solutions are simi- the CAPRI application program interface (API).
lar for both of the adaptation methods. The adaptive The initial coarse mesh shown in Fig. 10 is employed
procedure strongly clustered nodes at the leading and in the lift error prediction and grid adaptation proce-
trailing edges of the wing and lightly clustered nodes dure at 0.40 Mach an angle of attack of 2 deg. The
at the shock location on the upper surface. Feature- lift coefficient is plotted as a function of mesh size
based adaptations of this configuration in Ref. 7 and 8 in Fig. 11. The adaptive procedure has a user speci-
focused on the leading edge and shock locations, but fied error tolerance of 0.1 for the lift coefficient error.
not the trailing edge. The unifornfly refined lift calculation is shown with the
adapted grid lift calculation and error predictions on
Lift Adaptation - EET
the adapted grid in Fig. 11. The Richardson extrapo-
The EET (Energy Efficient Transport) 34 initial lation value is not shown because a reasonable linear
coarse mesh is shown in Fig. 10. The initial grid spac- fit of the last three points was not possible. The lift
ing distribution is specified manually to resolve the coefficient is calculated on an adapted mesh one-tenth
9 OF 11
AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2002 3286