Blonigan, Patrick

Carton De Wiart, Corentin

Diosady, Laslo T.

Ekelschot, Dirk

Garai, Anirban

Murman, Scott M.

NASA Technical Reports Server

RECENT DEVELOPMENTS FOR THE 
EDDY SOLVER
Scott M. Murman, Laslo T. Diosady, Anirban Garai, Corentin Carton 
de Wiart, Patrick Blonigan, Dirk Ekelschot
NASA Ames Research Center
https://ntrs.nasa.gov/search.jsp?R=20190027484 2019-09-26T19:15:23+00:00Z
Background
S. Murman 20-oct-15
• NASA CFD living off R&D investment 
from previous generation 
• Steady-state, complex geometry 
• Internal/external advocacy for new tools 
restarted fundamental R&D 
• Primarily numerical methods, turbulence 
modeling to-date 
• Use exascale computing to open new 
possibilities 
• Certification by simulation 
• Multi-disciplinary, multi-physics, robust error 
estimates, …
2
S. Murman 9-jan-18
Background
3
• Our group developing DG spectral-element capability for scale-
resolving simulations of separated flows over the past few years 
• Infrastructure known as the eddy solver 
• Provide overview of goals and technical approach 
• Fill in material not presented at AIAA
S. Murman 9-jan-18
• Complex geometry - unstructured mesh 
• Complex physics - scale-resolving methods 
• High-Re, combustion, chemistry - fully implicit methods 
• Computational intensive - high-order, adaptive methods 
• Multi-disciplinary, multi-physics - robust, extensible methods
Target Applications
4
OVERFLOW FUN3D
TURBODPLR
VTF
S. Murman 9-jun-16
Diagonalized ADI Preconditioner
5
A = (D+x ⌦ I ⌦A+x ) + (D x ⌦ I ⌦A x ) + (D+y ⌦ I ⌦A+y ) + (D y ⌦ I ⌦A y )
Upwind Jacobian
A = (I ⌦ I ⌦Rx)((D+x ⌦ I ⌦ ⇤+x ) + (D x ⌦ I ⌦ ⇤ x ))(I ⌦ I ⌦RTx )
+(I ⌦ I ⌦Ry)((D+y ⌦ I ⌦ ⇤+y ) + (D y ⌦ I ⌦ ⇤ y ))(I ⌦ I ⌦RTy )
Transform to characteristic variables and factor
• Can be simplified using mean speed over element 
• Approximate solve using 1D scalar problems in each direction 
• Convert to characteristic variables in X 
• Solve 5 1D scalar problems along lines in X 
• Repeat for Y, Z 
• Solve is dominant cost - use optimized matrix-matrix operations
• Matrix-free, Newton-Krylov solver, tensor-product
S. Murman 9-jan-18
Implicit Performance
• Same memory usage and better efficiency than explicit scheme  
• Implicit enables full space-time entropy-stable formulation 
• Current multi-physics work extending to monolithic framework
6
103 104 105 106 107
10−5
10−4
10−3
10−2
10−1
Work Units
Er
ro
r
 
 
N = 2
N = 4
N = 8
N = 16
0
1
2
3
4
5
0 1 2 3 4 5 6 7 8 9
C
os
t
re
la
tiv
e
to
ex
pl
ic
it
R
K
4
CFL = c th
N = 2
N = 4
N = 8
AIAA 2013-2870, JCP Vol. 330
S. Murman 9-jan-18
Turbulence Gradients
• Many groups following similar paths to develop scale-resolving 
capability 
• In order to meet goals need reliable/robust/general methods to 
calculate gradients 
• Adjoint methods for error estimation, adaptation, design, … 
• Turbulence is inherently chaotic and hard to model
7
S. Murman 9-jan-18
NACA 0012 Ma=0.2, α = 20o
8
0 10 20 30 40 50 60 70 80
Time
0.016
0.018
0.020
0.022
0.024
0.026
0.028
D
ra
g
0 20 40 60 80 100 120 140 160
Time
0.014
0.016
0.018
0.020
0.022
0.024
0.026
0.028
0.030
D
ra
g
0 20 40 60 80 100 120 140 160
Time
0.010
0.015
0.020
0.025
0.030
0.035
0.040
D
ra
g
260 280 300 320 340 360 380 400 420 440
Time
10 3
10 2
10 1
100
101
102
103
104
105
106
107
108
k k2
180 200 220 240 260 280 300 320 340 360
Time
0
50
100
150
200
250
300
k k2
50 60 70 80 90 100 110 120 130
Time
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
k k2
0.016 0.018 0.020 0.022 0.024 0.026 0.028
Drag(t)
0.016
0.018
0.020
0.022
0.024
0.026
0.028
D
ra
g(
t
+
 
t)
0.014 0.016 0.018 0.020 0.022 0.024 0.026 0.028 0.030
Drag(t)
0.014
0.016
0.018
0.020
0.022
0.024
0.026
0.028
0.030
D
ra
g(
t
+
 
t)
0.010 0.015 0.020 0.025 0.030 0.035 0.040
Drag(t)
0.010
0.015
0.020
0.025
0.030
0.035
0.040
D
ra
g(
t
+
 
t)
Re=800 Re=1600 Re=2400
Time
Phase
Adjoint log scalelinear scale
L. Diosady 03-01-17
T106 Low Pressure Turbine
9
• DNS - assume 
minimal modeling 
errors 
• Separation leading to 
transition/vortex 
shedding on suction-
side of blade 
• Fully turbulent wake
Re = 80,000, Min = 0.243, α = 32.7, Mout = 0.65 
L. Diosady 03-01-17
0 0.5 1 1.5
Time
1
1.1
1.2
1.3
1.4
1.5
Ax
ial
 F
or
ce
, = 32.7
, = 32.701
Sensitivity to Inflow Boundary Condition
10
• Modify inlet flow angle from α = 32.7 to α = 32.701
Domain flow-
through time
Convective disturbance
hits leading edgeAcoustic disturbance hits leading edge
S. Murman 9-jan-18
• Output is integrated axial 
force  
• Also define output without 
temporal normalization
Adjoint of mean Axial Force
11
J¯ =
1
T
Z T
0
Fx(u(⌧))d⌧
J(t) =
Z t
0
Fx(u(⌧))d⌧
Range: [-1e6, 1e6]
L. Diosady 03-01-17
Sensitivity computed using adjoint
12
0 0.5 1 1.5
Time
10 -8
10 -4
10 0
10 4
10 8
|J(
,
+"
 ,
)-J
(,
)|
Finite Difference
Adjoint
 J(t) = J(t;↵+ ↵)  J(t;↵) =
Z t
0
Fx(u(⌧ ;↵+ ↵))  Fx(u(⌧ ;↵))d⌧
⇡
Z t
0
 (⌧ ; t,↵)TR(u(⌧);↵+ ↵)
• Sensitivity computed 
using adjoint only valid 
for very short time 
windows 
• Adjoint computed using 
long time window blows 
up 
• Sensitivity computed 
using short time window, 
not representative long 
time behavior
13
Adjoint-based error indicator
S. Murman 20-oct-15
• Unbounded adjoint not useful for error estimation 
• Estimate is orders of magnitude larger than actual signal 
• Error localization simply flags regions where adjoint is large
• Estimate error using dual-
weighted residual method 
(Becker & Rannacher 1995) 
• Localize error 
✏ = J(u)  J(uH) ⇡ RH(uH , h)
✏ ⌘ RH(uH , h|)
S. Murman 9-jan-18
Adjoint growth with mesh resolution
14
0 0.5 1 1.5
Time
10 -1
10 1
10 3
10 5
10 7
10 9
|*
|
Coarse
Fine
• Refined mesh has essentially 
double mesh resolution near 
separation region 
• Increase mesh resolution 
results in faster growth of 
adjoint (i.e. larger Lyapunov 
exponent) 
• Adaptation mechanism is not 
convergent
S. Murman 9-jan-18
Minimum Turbulent Flow Unit
• Smallest channel that can sustain 
turbulent flow (Jimenez and Moin, 
1991).  
• Very good agreement with turbulent 
channel statistics below y+=40 
• Roughly 150 positive Lyapunov 
exponents 
• Compute sensitivity of total kinetic 
energy using adjoint shadowing methods
15
Q-Criterion	
isosurfaces	colored	
by	x-momentum
Blonigan et al. AIAA 2016-0296, JCP Vol. 348
S. Murman 9-jan-18
Least-Squares Shadowing Adjoint
• Shadowing adjoint does not exhibit exponential 
growth 
• Adjoint provides physical insights 
• Largest adjoint magnitudes occur before 
“bursting/blooming” of turbulence indicated by 
wall shear stress τ. 
16
t=12.13 t=13.18
Q-Criterion	isosurfaces	colored	by	x-momentum
Turbulence	“Bursting”
S. Murman 9-jan-18
Flow Unit Adjoint Field
• Integrated kinetic energy adjoint shows when and where flow is 
most susceptible to flow instabilities
17
Q-Criterion	isosurfaces	colored	
by	x-momentum
Adjoint	X-momentum	
isocontours	for	±2.0
Walls	colored	by	shear	
force	Magnitude
S. Murman 9-jan-18
• Adjoint	X-momentum	field	for	flow	unit	
prior	to	turbulence	“blooming”:
Optimal Perturbation for Transition
• Streamwise velocity magnitude contours 
for a flow perturbation optimized to 
increase the kinetic energy of Re=610 flow 
over a flat plate (Cherubini et al. 2010, 
JFM):
18
Contour	lines:	X-momentum	adjoint	
Color	map:	Z-vorticity
Solid	lines:	domain	length	=	400	units	
Dotted	lines:	domain	length	=	800	units
X-momentum	perturbations	suggested	by	the	adjoint	are	similar	
to	the	optimal	velocity	perturbations	computed	by	Cherubini	et	al.
S. Murman 9-jan-18
Scramjet Internal Flow
Adapted mesh with L1 norm 22,566 vertices:
39 Anisotropic Mesh Adaptation
Figure 17: Variational Hex-Dominant Meshing and comparison with [Mare´chal 2009] (data kindly provided by paper’s author).
5.4 Variational Hex-Dominant Meshing
Our framework extends naturally to
variational hex-dominant meshing. The
volume framework is used (i.e., Lp-
CVT is computed on the clipped
Voronoi diagram), with an anisotropy
matrix M associated with each inte-
gration simplex. The anisotropy is
attached to the facets of the surface,
defined as in the previous section,
and each integration simplex uses the
anisotropy Mf of the facet f nearest
to its centroid, queried using Approx-
imated Nearest Neighbor [Mount and
Arya 1997]. Note that the Delaunay tri-
angulation restricted to the interior of
S is smaller than the interior of S, as
shown in the opposite figure. To avoid
this “shrinkage”, vertices on the bound-
ary require a special treatment, simi-
larly to [Tournois et al. 2009]. Vertices on the boundary are com-
puted by the surface framework, as in the previous subsection, us-
ing the same values ofMf . These vertices are characterized by a
Figure 18: Mesh repair, comparison with Delaunay refinement
[Busaryev et al. 2009], (data kindly provided by paper’s author).
non-empty intersection of their Voronoi cells with the surface S.
After the optimization step, the algorithm computes the Delaunay
triangulation of the vertices restricted to the interior of S, using the
combinatorial information computed in Section 4.1. Each vertex of
configuration D (Voronoi vertex inside S) corresponds to a tetrahe-
dron. Finally, the tetrahedra are merged to form hexes, using com-
binatorial optimization [Meshkat and Talmor 2000]. The subgraphs
of the simplex graph that correspond to hexes are pattern-matched
and hexes are generated, in priority order, determined by dihedral
angles and face planarity.
Figure 17 shows a hex-dominant mesh (boundary and cross-
sections) of a CAD model with timing and stats, and compares the
result with an octree-based method [Mare´chal 2009]. The advan-
tage of the octree-based method are that it generates a pure hex-
ahedral mesh with a nearly regular pattern and that it is fast. The
advantages ofLp-CVT is that boundaries are well respected regard-
less their orientation, and that it generates homogeneous element
sizes. Figures 1 and 19 show examples with smooth surfaces.
5.5 Discussion - Limitations
Lp-CVT is robust to T-junctions and artificial borders (see the em-
pirical results in the paper and supplemental material), i.e. meshes
that have a correct geometry but a non-conforming discretization.
However, Lp-CVT might not handle properly models that are ge-
ometrically inconsistent. Figure 18 compares the result with the
CADmodel repair algorithm proposed in [Busaryev et al. 2009] and
illustrates this limitation of Lp-CVT. The input data has some gaps.
Lp-CVT misses triangles when a Voronoi edge passes through a
gap (more triangles are missed when the number of vertices is in-
creased). In practice, hole filling heuristic fixes the problem in most
cases. More formally, we think that the "-sampling formalism and
its extensions [Boissonnat and Oudot 2006] provide a mean of de-
veloping the theoretical tools to study how to recover the missing
triangles with theoretical guarantees. Note also that the continuity
of FLp remains to be studied. Finally, we also mention that under
extreme anisotropy (100x ratio between the lengths of axes) con-
vergence becomes much slower because of the conditioning of the
Hessian of FLp .
119:8       •       B. Lévy et al.
ACM Transactions on Graphics, Vol. 29, No. 4, Article 119, Publication date: July 2010.
Metric-based Mesh Adaptation
• Meshing/adaptation uses similar strategy as solver development 
• Design from scratch to meet objectives 
• Automatic, hex-dominant, feature-aligned, … 
• Based on mathematically robust / provable algorithms 
• Align using Riemannian metric field 
• Error estimates, surface curvature, …  
• Reduces dimension of problem, e.g. 3D to 2D
19
Scramjet Loseille et al.
Lévy & Liu
S. Murman 9-jan-18
Bounded Lp-CVT Mesh Adaptation
• Extend Lévy & Liu approach to bounded/finite domain 
• Hierarchical approach: edges -> surfaces -> volumes 
• Preliminary 2D proof of concept
20
Ekelschot et al. AIAA 2018-1501
-2 -1 0 1 2
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
S. Murman 9-jan-18
(a) Convergence of the energy value. (b) Convergence of the norm of the gradient.
Fig. 9 Convergence history for the metric given in Eqn. (13)
computational domain in this case is ⌦ = [ 1, 1]2 and the function given in Eqn. (14) is plotted in Fig. 10. We start
Fig. 10 Metric resembling a weak s-shock.
by measuring the length of the boundary edges subject to the underlying metric and we determine the number of
boundary nodes that are required to mesh each boundary edge such that each sub-edge has unit length. The location of
these boundary vertices are then found using the Newton method. Once, the boundary edges are meshed in a metric
conforming way, we randomly place 150 internal vertices and use that as an initial guess for our optimization procedure.
The initial state with its Delaunay triangulation is plotted in Fig. 11a and the corresponding Voronoi diagram is given in
Fig. 11b. The minimal energy state that is obtained after running the proposed Lp-CVT algorithm for the given initial
point distribution and the given underlying metric is presented in Fig 11c and 11d. The convergence history is shown in
16
Bounded Lp-CVT Mesh Adaptation
• Extend Lévy & Liu approach to bounded/finite domain 
• Hierarchical approach: edges -> surfaces -> volumes 
• P eliminary 2D proof of concept
21
Ekelschot et al. AIAA 2018-1501
S. Murman 9-jan-18
Multi-physics Approach
• Each module is a separate 
physics/set of equations to solve 
• Physics are coupled 
• General approach is required 
• Shock capturing 
• Chemistry 
• Combustion 
• …
22
Wall-model
Cable	solver
Fabric	solver
Non-reflecting	BC
Deforming	mesh	solver
Wake	LES	model
6-DOF
S. Murman 9-jan-18
Multi-physics Capability
23
• General monolithic implicit multi-physics solver 
• Exploits software-design to modularize physics, 
discretization, etc. without sacrificing efficiency 
• Similar to partitioned approach but w/o loss of conservation, 
accuracy, stability 
• Automatically support primal, adjoint, and tangent equations 
for any system 
• All leverage same optimized kernels, solver, etc. 
• Supports CG, DG, C1-DG discretizations 
• Easily extends, e.g. HDG, optimized basis for b.c. 
• Does not require researcher to understand entire code to 
leverage 
Carton de Wiart et al. AIAA 2018-1400
S. Murman 9-jan-18
Current Status
• Synergy between R&D led by RCA and engineering projects 
• Internal focused engineering partnerships 
• Turbomachinery (ARMD AATT) 
• Transonic buffet (SLS and NESC) 
• Parchute FSI (Orion CPAS, STMD, and NESC) 
• Aft-body aeroheating and JI (Orion CPAS, STMD) 
• External collaborations 
• Public domain license 
• MIT, Michigan, Stanford, UIUC, UNM, U. Colorado, UTIAS, 
Boeing, Cenaero 
• Currently supporting 4 PhD projects
24
S. Murman 9-jan-18
Backup
25
S. Murman 9-jan-18
Context
• NASA has healthy infrastructure of CFD tools 
• OVERFLOW, FUN3D, Cart3D, Loci/CHEM, DPLR 
• Primarily FD/FV, RANS-based technology 
• Goal is to complement existing suite, not replace 
• Unsteady complex physics, e.g. separation, shock/BL interaction, … 
• Use mathematically robust algorithms and procedures 
• Required for error estimates, uncertainty quantification, and 
automation 
• Achieve efficiency by improved methods, not short-cuts
26


Recent Developments for the Eddy Solver

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server