MARCH/APRIL 2007 THIS ARTICLE HAS BEEN PEER-REVIEWED. 21
Stochastic Computational
Fluid Mechanics
S TOCHASTIC
M ODELING
N
umerical solution of differential
equations with stochastic inputs has
direct impact on simulations of flow
systems in at least three areas: uncer-
tainty quantification, stability of noisy systems, and
coarse-grained and multiscale representation. In
the first category, the past few years have seen an
increasing interest in uncertainty quantification in
large-scale numerical simulations. In simulations,
just as in experiments, we often question our
results’ accuracy and construct error bounds a
posteriori, but the new objective is to model un-
certainty from the simulation’s start, not to do it
simply as an afterthought. Researchers have used
numerical accuracy and error control in flow sim-
ulations for some time now, at least in the more
modern discretizations, but there’s still an uncer-
tainty component associated with the physical
problem (specifically, with diverse factors such as
constitutive laws, boundary and initial conditions,
transport coefficients, source and interaction terms,
geometric irregularities, and so on).
We encounter noisy nonlinear systems in many
applications, from the nanoscale (as in self-assem-
bly processes) to the macroscale (such as large, sud-
den, disturbances in flow past an aircraft).
Bifurcations and chaotic transitions in stochastic
dynamical systems can differ greatly from instabil-
ities in deterministic dynamical systems.
1
How ex-
actly extrinsic stochasticity interacts with intrinsic
stochasticity caused by the systems’ nonlinear in-
teractions is an intriguing and not-well-understood
problem. In turbulent boundary layers, in which a
wide range of small scales exists, the mean flow
seems to be totally unaffected, even in the presence
of substantial background turbulence.
2
In contrast,
for other flow systems at a low Reynolds number,
even small amounts of noise can have a significant
effect on the mean flow’s structure.
3
The third area in the aforementioned list in-
volves complex systems with an extremely large
number of degrees of freedom—models of turbu-
lent flows at a very high Reynolds number or atom-
istic simulations of mesoscopic size systems.
Frequently, we use a coarse-grained procedure to
remove degrees of freedom making relatively small
contributions to the overall system’s energy. The
removed degrees of freedom are usually accounted
for by stochastic terms in the dimensionally re-
duced “effective” or “upscaled” equations. In the
dissipative particle dynamics method,
4
for example,
coarse-graining of the molecular dynamics method
Stochastic simulations in computational fluid dynamics let researchers model uncertainties
beyond numerical discretization errors. The authors present examples of stochastic
simulations of compressible and incompressible flows and provide analytical solutions for
verifying these newly emerging methods for stochastic modeling.
GUANG LIN, XIAOLIANG WAN, CHAU-HSING SU,
AND GEORGE EM KARNIADAKIS
Brown University
1521-9615/07/$25.00 © 2007 IEEE
Copublished by the IEEE CS and the AIP
22 COMPUTING IN SCIENCE & ENGINEERING
leads to a system of stochastic ordinary differential
equations (ODEs) that must be solved efficiently
for a very large number of particles.
In this article, we present an efficient nonstatis-
tical approach for modeling extrinsic stochasticity
for both compressible and incompressible flows.
We consider two prototype cases—flow past a
wedge and flow past a bluff body—and quantify un-
certainties associated with inflow disturbances, ran-
dom motions, and random geometric roughness.
This is a new field, so properly verifying the results
is particularly important. To this end, we’ve derived
analytical solutions for small disturbances based on
stochastic perturbation analysis, and we employ
these solutions in the verification procedure. Vali-
dation of stochastic computational fluid dynamics
(CFD) models is a much more difficult task because
it requires very accurate characterization experi-
mentally of both inputs and outputs, which is cur-
rently lacking in existing data.
Polynomial Chaos and Its Variants
Polynomial chaos represents a stochastic process
by a spectral expansion based on Hermite or-
thogonal polynomials in terms of Gaussian ran-
dom variables. Roger G. Ghanem and Pol D.
Spanos pioneered its use in solving stochastic dif-
ferential equations by employing a Galerkin pro-
jection to derive an equivalent system of
deterministic equations; this can, typically, be
solved with standard numerical techniques.
5
Dongbin Xiu and George E. Karniadakis devel-
oped generalized polynomial chaos (gPC), which
uses a broader family of trial bases and is based on
the orthogonal polynomials from the Askey
scheme.
6
As an example, gPC can express a gen-
eral second-order random process T() as
, (1)
where is the random event, and the family {
i
} is
an orthogonal basis in terms of the random vector
() with the following orthogonality relation
,
where
ij
is the Kronecker delta, and denotes the
ensemble average with respect to the probability
density function (PDF) of . For a certain random
vector , we can choose the gPC basis {
i
} in such a
way that its weight function has the same form as
the PDF of . Gaussian random variables, for ex-
ample, are associated with Hermite polynomials,
uniform random variables are associated with
Legendre polynomials, and so forth.
Global gPC expansions work effectively for
many physical applications, but increasing the
polynomial order might not be efficient for some
problems, such as those with random frequencies
or low regularity in the parametric space. Based on
gPC, Xiaoliang Wan and Karniadakis developed
the multi-element extension (ME-gPC), which de-
composes the random space into finite elements as
in the deterministic finite element method.
7
Fig-
ure 1 illustrates this idea in the one-dimensional
case. We assume that the range [a, b] of the 1D
random variable is decomposed into elements e
k
:= [a
k
, b
k
]. We then define a new random variable
k
, k = 1, 2, …, N, in each random element, e
k
, as
, (2)
with a re-scaled PDF
, k = 1, …, N,
where N is the number of random elements, f() is
the PDF of , and
is the probability that is located in random ele-
ment e
k
. We first approximate the desired random
field u() locally via gPC, where the degree of per-
turbation is effectively decreased by Equation 2
from O(1) to O((b
k
– a
k
)/2)). Subsequently, we
gather the information from all random elements
to obtain the statistics of u().
Because the PDF of is also decomposed to-
gether with the random space, gPC’s orthogonal-
ity in the entire random space will be, in general,
lost in random elements. However, given an arbi-
trary PDF, we can construct the following orthog-
Pr( ) ( )
ξξξ
∈=
∫
efd
k
e
k
f
f
e
ba
kk
k
k
kk
()
(( ))
Pr( )
ζ
ξζ
ξ
=
∈
−
2
ξζ
=
−
+
+ba ba
kk
k
kk
22
φφ φ δ
ij i ij
, =
2
TT
ii
i
()
ˆ
(( ))
ωφξω
=
=
∞
∑
0
f()
aa
k
b
k
b –1 1
e
k
e
k
f
k
(
k
)
Figure 1. Schematic of decomposition and mapping in multi-
element generalized polynomial chaos (ME-gPC). The procedure is
similar to the deterministic finite element method.
MARCH/APRIL 2007 23
onal system numerically on the fly, from the three-
term recurrence relationship
i+1
(
k
)= (
k
–
i
)
i
(
k
) –
i
i–1
(
k
), i = 0, 1, ,
0
(
k
) = 1,
–1
(
k
) = 0,
where {
i
(
k
)} is a set of (monic) orthogonal poly-
nomials,
i
(
k
) = + lower-degree terms, i = 0, 1, … ,
and the coefficients
i
and
i
are uniquely deter-
mined by a positive (probability) measure f
k
(
k
).
The orthogonal system {
i
(
k
)} will then serve as a
local gPC basis. For problems related to low regu-
larity in the parametric space, we use the ME-gPC
method adaptively to achieve high efficiency.
8
Galerkin versus Collocation Projection
For simplicity and clarity, we consider the stochas-
tic equation L(x, t, (); u) = f (x, t; ()) with a gen-
eral (nonlinear) differential operator L, where x
d
, d = 1, 2, 3, indicates the physical space and t the
time. In the Galerkin formulation, we first apply
the gPC expansion to
and
.
Here P is the total number of basis modes. Then,
performing the Galerkin projection on both sides
of the equation, we obtain
, (3)
where k = 0, 1, …, P – 1. In contrast, in the collo-
cation formulation, we employ Delta functions
( –
k
) as test functions, k = 0, ..., M – 1, where
{
k
} is a proper set of grid (quadrature) points on
the range of (), and M is the number of grid
points. By applying collocation projection on both
sides of the equation, we obtain
L(x, t,
k
; u) = f(x, t;
k
). (4)
Obviously, Equation 4 shares the same format as
the original equation, whereas Equation 3 is dif-
ferent. In particular, if the operator L is nonlinear,
Equation 3 is a system of ordinary or partial differ-
ential equations with the unknowns u
i
coupled to-
gether, which introduces a real challenge to design
an efficient deterministic solver. For complex fluid
dynamical systems, such as compressible flows, it’s
easier to use the collocation projection to obtain
the governing equations.
Note that the stochastic collocation method is
equivalent to solving a deterministic problem at
each grid point. Once we obtain the numerical so-
lutions at all collocation points, we can evaluate the
random solution’s statistics using the correspond-
ing quadrature rule, for example,
,
where w
k
represents the integration weights, and
denotes the expectation.
High Dimensionality and Sparse Grids
Xiu and Jan S. Hesthaven
9
constructed a stochas-
tic collocation method based on sparse grids by
using the Smolyak algorithm.
10
This algorithm is
a linear combination of tensor product formulas,
and the resulting grid set has a significantly
smaller number of grids compared to the full ten-
sor product rule. Sparse grids depend weakly on
the random space’s dimensionality, and hence are
more suitable for applications with large dimen-
sional random inputs. Recently, sparse grids for
high-dimensional integrals have also received a lot
of attention.
11–13
Representation of Stochastic Inputs
We represent the stochastic inputs with the
Karhunen-Loeve (KL) decomposition
14
given by
, (5)
where v(x; ) denotes the random process,
denotes the corresponding mean, {
i
()} is a set of
uncorrelated random variables with zero mean and
unit variance, and x is the spatial or temporal co-
ordinate. Also,
i
(x) and
i
are the eigenfunctions
and eigenvalues of the covariance kernel R
hh
(x, y),
respectively—that is,
D
R
hh
(x, y)
i
(y)dy =
i
i
(x). (6)
Specifically, for a time-dependent stochastic
process, we usually consider the exponential co-
variance kernel
v()x
vv
i
i
ii
(; ) () () ( )xx x
ωλψξω
=+
=
∞
∑
1
[ ]( ,) ( ,, ) ( ) ( ,, )u t ut f d ut w
kk
k
M
xx x=≈
=
−
∑
ξξξ ξ
0
1
∫∫
〈
⎛
⎝
⎜
⎞
⎠
⎟
〉=〈 〉
=
−
∑
Lt u f
ii
i
P
kk
x,, ( ); , ,
ξω φ φ φ
0
1
ft f
i
P
ii
(,;( ))x
ξω φ
=
=
−
Σ
0
1
ut u
i
P
ii
(,;( ))x
ξω φ
=
=
−
Σ
0
1
24 COMPUTING IN SCIENCE & ENGINEERING
, (7)
where A is the correlation length. This corresponds
to a first-order Markov chain in time, which in dis-
crete form is
v
0
=
0
, v
i+1
= cv
i
+ q
i+1
, i 0, (8)
where , , and t is the time step.
Discrete forms of stationary spatial-dependent
random processes usually differ from time-depen-
dent ones. This is due to the time series’ unilateral
nature, which only depends on past values as op-
posed to the spatial process’s dependence in all di-
rections. By solving the following stochastic
Helmholtz equation, we can obtain the correspond-
ing spatial covariance kernel R
hh
(v(x
1
, ), v(x
2
, )),
15
v – k
2
v = f(x), (9)
where the random forcing term f(x) is a white-
noise process, a function of the spatial vector x
that satisfies
[f(x
1
)f(x
2
)] = (x
1
– x
2
). (10)
Equation 9 corresponds to a discrete model for
a second-order autoregressive process in space,
where the value of v at x depends on its neighbors.
Here’s a simple 1D example:
. (11)
The parameter c is the constant correlation coeffi-
cient, a gives the measure of the stochastic forcing’s
strength, and
i
are independent random variables
with zero mean and unit variance.
Stochastic Simulations
Let’s first look at some results for compressible in-
viscid steady flows and then for incompressible un-
steady flows.
Stochastic Compressible Flows
For supersonic flows, random inflow disturbances
can substantially modify the flow field, thus ren-
dering traditional predictive CFD tools invalid.
We study the shock dynamics in two-dimensional
supersonic flows past a wedge, assuming two dif-
ferent random disturbances. These correspond to
random inflow velocity or random oscillations of
the wedge around its apex, and they can be steady
in time or time-dependent (the latter could be a
model for aeroelastic motions, such as flutter os-
cillations). Figure 2 shows a schematic. We con-
sider small perturbations and assume that the
shock slope’s perturbation is small. We approxi-
mate the flow between the shock and the wedge
as isentropic.
We denote the wedge angle by
0
, the shock an-
gle by
0
, and the incoming flow velocity W
1
with
its normal component u
1
= W
1
sin
0
. The stream-
lines behind the shock are parallel to the wedge
surface, and we denote the velocity by W
2
and its
normal component to the shock by u
2
= W
2
sin(
0
–
0
) = W
1
cos
0
tan(
0
–
0
).
First, we assume that the wedge angle is fixed,
but the shock is subject to random inflow pertur-
bation that we describe as a uniform random vari-
able. We’ve obtained an analytical solution for the
perturbed shock path
, (12)
where s = tan(
0
–
0
), M
1
is the perturbed part
of the inflow Mach number, and H(
0
) = dM
1
/d.
The mean and variance of the perturbed shock
path are then
z(x; ) = 0,
zx x s
M
H
(;) ( )
()
ξ
χ
=+
Δ
1
2
1
0
v
c
vv a
iii i
=++
+−
2
11
()
ξ
qc=−1
2
ce
t
A
=
−Δ
Rvt vt
vt vt e
hh
((,),(,))
[( , )( , )]
12
12
ωω
ωω
==
−−−||tt
A
12
y
x
2
p
2
c
2
1
p
1
c
1
W
1
W
2
→
→
0
u
1
u
2
0
0
v
u
v
2
v
1
Random perturbed shock
Unperturbed shock
M
1
> 1
d
0
–
0
Figure 2. Sketch of supersonic flow past a wedge and definition of a
coordinate system. We assume two different random disturbances
corresponding to random inflow velocity or random oscillations of
the wedge around its apex. We’re interested in investigating how
such random disturbances affect the shock path.
MARCH/APRIL 2007 25
(13)
Next, we assume that the wedge inflow is deter-
ministic and describe the wedge oscillations as a
uniform random variable. In this case, we obtain
. (14)
The mean and variance for this case are similar to
the previous case, consistent with physical intuition.
Numerically, we solve the 2D Euler equations
following the polynomial chaos approach for ran-
dom inflow and random oscillations. In the latter
case, we use a transformation based on a boundary-
fitted coordinate system approach, so that we solve
the Euler equations in a stationary domain:
u
t
+ f(u)
x
+ g(u)
y
= T(u), (15)
where u = [, m, n, E]
T
, f = [u, um + p, un, u(p + E)]
T
,
and g = [v, vm, vn + p, v(p + E)]
T
; for random inflow,
T = 0, whereas for random oscillations, T = [0,
–cos
0
(u
w
)
t
, sin
0
(u
w
)
t
, (u
w
)
t
(nsin
0
– mcos
0
)]
T
.
Here, u
w
represents the wedge’s random motion.
In these examples, we consider the case in which
the inflow velocity is perturbed by random fluctu-
ation, which we describe as a random variable with
amplitude = 0.01 and = 0.18. Figure 3 presents
the mean of density (x, y; ) corresponding to in-
flow perturbation described as a random variable
with amplitude = 0.18. We note that the mean re-
sembles a fan expansion. Figure 4 plots the per-
turbed shock path’s variance as a function of the
distance from the wedge apex x on the wedge sur-
face for two amplitude values: = 0.01 and = 0.18.
From Equation 13, we know that the perturbed
shock’s variance is proportional to x
2
. Indeed, for
amplitude = 0.01, we can verify that the pertur-
bation solution from Equation 13 and the numer-
ical solution for amplitude = 0.01 match exactly.
However, for amplitude = 0.18, we see in Figure
4 that the numerical solution deviates from the
analytical solution in Equation 13 because that
equation doesn’t hold for large amplitudes.
Stochastic Wedge Roughness
Next, let’s consider the perturbation of an oblique
attached shock in supersonic flow past a wedge due
to random roughness on the wedge’s surface.
16
Fig-
ure 2 defines the problem, and we use the same as-
zx x s
M
H
(;) ( )
()
ξ
εξ
χ
=− +1
2
1
0
Var z x x s
M
H
xs
((;)) ( )
()
(
ξ
εξ
χ
=+
〈〉
=+
222
1
22 2
2
0
2
1
1
222
1
22
2
0
3
)
()
.
M
H
ε
χ
Deterministic
2.2
2.1
2.0
1.9
1.8
1.7
1.6
Figure 3. Mean of (x, y; ) induced by random inflow perturbation
described as a random variable with amplitude, = 0.18. Note that
the mean resembles a fan expansion.
x
<Z
2
>/
2
012345
0
5
10
15
Perturbation solution
ME-gPC: = 0.01
ME-gPC: = 0.18
Figure 4. Variance of perturbed shock paths as a
function of x induced by random inflow
perturbation, which is described as a random
variable with amplitude = 0.01 and = 0.18.
The perturbed shock path’s variance obtained
from the numerical simulations for amplitude
= 0.01 is proportional to x
2
, which matches with
the perturbation solution. However, for
amplitude = 0.18, the numerical solution
deviates from the analytical solution because
the perturbation analysis doesn’t hold for large
amplitude.
![Figure 7. Schematic of the domain and mesh for noisy flow past a circular cylinder. The size of the domain is [–15D, 25D] [–9D, 9D], and the cylinder is at the origin with diameter D = 1. The mesh](/figures/figure-7-schematic-of-the-domain-and-mesh-for-noisy-flow-1j8uy1ed.webp)
