A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media

1
A stabilized stochastic finite element
second-order projection method for modeling
natural convection in random porous media
Xiang Ma and Nicholas Zabaras
Materials Process Design and Control
Laboratory Sibley School of Mechanical and
Aerospace Engineering188 Frank H. T. Rhodes
Hall Cornell University Ithaca, NY
14853-3801 Email zabaras_at_cornell.edu URL
http//mpdc.mae.cornell.edu
Materials Process Design and Control Laboratory
2
Outline of the presentation

• Motivation and problem definition
• A stabilized second-order projection method
• Representation of stochastic processes
• Stochastic finite element method formulation
• Numerical examples

Materials Process Design and Control Laboratory
3
Transport in heterogeneous media
- Thermal and fluid transport in heterogeneous
media are ubiquitous, e.g. solidification -
Range from large scale systems (geothermal
systems) to the small scale - Complex phenomena -
How to represent complex structures? - How to
make them tractable? - Are simulations
believable? - How does uncertainty propagate
through them?

• To apply physical processes on these
  heterogeneous systems
• worst case scenarios
• variations on physical properties

Materials Process Design and Control Laboratory
4
Problem of interest
0, u v 0
Natural convection in random porous media
1 u v 0
0 u v 0
Deterministic governing equation
0, u v 0
Pr -Prandtl number Ra- Raleigh number Da- Darcy
number
Materials Process Design and Control Laboratory
5
Problem of interest
0, u v 0
Natural convection in random porous media
1 u v 0
0 u v 0
Stochastic governing equation
0, u v 0
Modeling porosity as a Gaussian random
field
Materials Process Design and Control Laboratory
6
Problem of interest Some difficulties
This model is important in alloy solidification
process. To understand the uncertainty
propagation in this process can help understand
error in the experiment results.
very hard to extend to stochastic formulation due
to the high degree of freedom. Especially, when
dealing with real porous media, it results in a
very high stochastic dimension.
A SUPG/PSPG formulation has been proposed to
solve this complex problem. ( Deep Zabaras,
2004). However, due to its coupling between
velocity and pressure, it resulted in a very
large linear system.
Materials Process Design and Control Laboratory
7
Problem of interest Some difficulties and
solutions

• The projection method for the incompressible
  Navier-Stokes equations, also known as the
  fractional step method or operator splitting
  method, has attracted widespread popularity.
• The reason for this lies on the decoupling of
  the velocity and pressure computation.
• It was first introduced by Chorin, as the
  first-order projection method ( also called
  non-incremental pressure-correction method).
  Later, it was extended to the second-order scheme
  ( also called incremental pressure-correction
  scheme ) in which part of the pressure gradient
  is kept in the momentum equation.
• Le Maitre et.al(2002) have developed a
  stochastic projection method to model natural
  convection in a closed cavity based on first
  order method.
• It is really helpful when solving the high
  stochastic dimension problem.

Materials Process Design and Control Laboratory
8
Problem of interest Some difficulties and
solutions

• Let us review first order method is

solve for intermediate velocity

projection step
eliminate
Once a finite element discretization is
performed, the matrix form is
eliminate

• The term may be understood
  as the difference between two discrete Laplacian
  operators computed in difference manner. This
  matrix turns out to be positive semi-definite,
  which increases the stability of the numerical
  method, thus allows equal-order interpolation.

Materials Process Design and Control Laboratory
9
Problem of interest Some difficulties and
solutions

• For the porous media problem considered here, due
  to the porosity dependence of the pressure
  gradient term in the momentum equation, we cannot
  omit the pressure term as is the case in the
  first-order projection method. We need to use
  second-order scheme.

porosity dependence

• However, if using FEM, second-order projection
  method is known that it exhibits spatial
  oscillation in the pressure field if not using a
  mixed finite element formulation for velocity and
  pressure. (J.L.Guermond et al. 1998)

Materials Process Design and Control Laboratory
10
Problem of interest Some difficulties and
solutions
pressure field in a lid-driven cavity problem
with stabilization (left) and without
stabilization (right).

• In order to utilize the advantage of the
  incremental projection method which retains the
  optimal space approximation property of the
  finite element and allows equal-order finite
  element interpolation, Codina and co-workers
  developed a pressure stabilized finite element
  second-order projection formulation for the
  incompressible Navier-Stokes equation. (Codina et
  al.2000)

• The method consists in adding to the
  incompressibility equation the divergence of the
  difference between the pressure gradient and its
  projection onto the velocity space, both
  multiplied by algorithmic parameters defined
  element-wise.

Materials Process Design and Control Laboratory
11
Pressure stabilized formulation

• We take these parameters as

where is the viscosity, is the local size
of the element , is the local velocity
and and are algorithmic constants, that
we take as and for linear
elements.

• Accordingly, the continuity equation is modified
  as follows

where is the stabilized parameter as
discussed before and the new auxiliary variable
is the projection of the pressure gradient
onto the velocity space.
Materials Process Design and Control Laboratory
12
Pressure stabilized formulation
Lets see how it works
weak form
A matrix version of this form is
Eliminating from equation, it is found
that
which is similar to

• So this stabilized method mimics the stable
  effect of the first-order method. Thus it allows
  the equal-order interpolation and stabilize the
  pressure.

Materials Process Design and Control Laboratory
13
Stabilized second-order projection method

• Step 1 Solve for the intermediate velocity
  in the momentum equation

• Step 2 Projection step

• Step 2 update step

Fully decoupled, can be solved one by one.
Materials Process Design and Control Laboratory
14
A numerical example

• The computation region is a square
  domain.
• The dimensionless length is taken as
• The wall porosity is taken as 0.4. The
  porosity increases linearly from 0.4 at the wall
  to 1.0 (pure liquid) at the core.
• The Rayleigh number is , Prandtl number
  is 1.0 and the Darcy number is
• The problem was solved with a mesh discretization
  
• and the time step was chosen as
• The simulation was run up to non-dimensional time
  1.0.

1 u v 0
0 u v 0
Free fluid
Porous Medium
Materials Process Design and Control Laboratory
15
A numerical example
First-order solution
Second-order solution
The Streamline pattern is symmetric
Isothermal
Streamline
Materials Process Design and Control Laboratory
16
A numerical example
First-order solution
Pressure
v velocity
u velocity
Second-order solution
The velocities are mainly distributed in the free
fluid region
Materials Process Design and Control Laboratory
17
Representation of stochastic processes

• For special kinds of stochastic processes that
  have finite variance-covariance function, we have
  mean-square convergent expansions

Series expansions

• Known covariance function

• Unknown covariance function

Karhunen-Loeve
Generalized polynomial chaos

• Best approximation in mean-square sense
• Useful typically for input uncertainty modeling

• Can yield exponentially convergent expansions
• Used typically for output uncertainty modeling

Materials Process Design and Control Laboratory
18
Karhunen-Loeve expansion

• is the mean of the process and
  let denote its covariance
  function, which is needed to be known a priori.
  
• is a set of i.i.d. random variables,
  whose distribution depends on the type of
  stochastic process.
• and are the eigenvalues and
  eigenvectors of the covariance kernel. That is
  ,they are the solution of the integral equation
• In practice, we truncate KLE to first M terms.

Materials Process Design and Control Laboratory
19
Generalized polynomial chaos

• Generalized polynomial chaos expansion is used
  to represent the stochastic output in terms of
  the input

Stochastic input
Askey polynomials in input
Stochastic output
Deterministic functions

• Askey polynomials type of input stochastic
  process
• Usually, Hermite (Gaussian), Legendre (Uniform),
  Jacobi etc.

Materials Process Design and Control Laboratory
20
Generalized polynomial chaos

• The polynomial chaos forms a complete orthogonal
  basis in the of the Gaussian random
  variables, i.e.

where denotes the expectation or ensemble
average. It is the inner product in the Hilbert
space of the Gaussian random variables
where the weighting function has the
form of the multi-dimensional independent
Gaussian probability distribution with unit
variance.
1 dimension
2 dimensions
n dimensions
Materials Process Design and Control Laboratory
21
PC computation
PC computation (B. Debusschere et al. 2004, Ma
Zabaras, 2007)

• Quadratic product

• Consider two random variables, and , with
  their respective GPCE

• We want to find the coefficients of

• The coefficient are obtained with a
  Galerkin projection, which minimizes the error
  of the resulting PC representation within the
  space spanned by the basis function up to order
  P.

Materials Process Design and Control Laboratory
22
PC computation

• Quadratic product

• Thus we can obtain

with

• Triplet product

• A similar procedure could also be used to
  determine the product of
• three stochastic variables

with
Materials Process Design and Control Laboratory
23
PC computation
Sampling approach for general nonpolynomial
function evaluations

• We consider a general nonlinear function
  , where w is

We need to express this function as

• Using the same method as before (Galerkin
  projection), we have

High dimension integration, use Latin Hypercube
sampling.
Thus we obtain
Materials Process Design and Control Laboratory
24
Stochastic Finite Element Formulation

• Since there are nonlinear functions of
  in the governing equation, we need to first
  express them in the polynomial basis using the
  method discussed before

• Since the input uncertainty is taken as a
  Gaussian random field, we use Hermite polynomials
  to represent the solution

Materials Process Design and Control Laboratory
25
Stochastic Finite Element Formulation

• Substitute the above equation into Eqs
  (1)-(5)and then performing a Galerkin projection
  of each equation by

Materials Process Design and Control Laboratory
26
Stochastic Finite Element Formulation

• In the non-linear drag term
  
  
• we assume that the magnitude of the velocity
  is determined by the
• mean velocity

• From the pressure update equation
• the stabilized parameter is determined by the
  kth coefficient of the spectral expansion. So we
  denote it as to emphasize that it is a
  function of the k-th coefficient

Materials Process Design and Control Laboratory
27
Stochastic Finite Element Formulation

• It is time consuming to evaluate the fourth
  order product term
• 
  
• directly using the method discussed before. To
  simplify this calculation, we introduce an
  auxiliary random variable as follows
• such that
• So our term of interest can be reduced to
• This form is now easier to calculate, since it
  only evolves third-order product terms, which is
  pre-computed.

Materials Process Design and Control Laboratory
28
Some computational details

• All calculations are performed using numerical
  algorithms provided in PETSc. (Portable,
  Extensible Toolkit for Scientific computation).
  The default Krylov solver was used for the linear
  solver).
• For the numerical computation of KLE, the SLEPc
  is used. ( Scalable Library for Eigenvalue
  Problem Computation). It is based on the PETSc
  data structure and it is highly parallel.
• In order to further reduce the computational
  cost, we solve each expansion coefficient
  individually. That means we split the systems of
  equations into (P1)(3d2) deterministic scalar
  equations.
• The computation is performed using 30 nodes on a
  local Linux cluster.

Materials Process Design and Control Laboratory
29
Natural convection in heterogeneous media large
dimension
Alloy solidification, thermal insulation,
petroleum prospecting Look at natural convection
through a realistic sample of heterogeneous
material
0, u v 0
1 u v 0
0 u v 0
0, u v 0

• The porosity is modeled as a Gaussian random
  filed. And the correlation function is extracted
  from a real porous media sandstone 1. The mean
  is 0.6.

1. Reconstruction of random media using Monte
Carlo methods, Manwat and Hilfer, Physical Review
E. 59 (1999)
Materials Process Design and Control Laboratory
30
Natural convection in heterogeneous media
Material Sandstone
Experimental correlation for the porosity of the
sandstone. Eigen spectrum is peaked. Requires
large dimensions to accurately represent the
stochastic space KLE is truncated after first 6
terms. So the stochastic dimension is 6.
Correlation function
Numerically computed Eigen spectrum
Materials Process Design and Control Laboratory
31
Natural convection in heterogeneous media
Eigenmodes for a 2-D domain
Mode 2
Mode 1
Mode 4
Mode 6
Materials Process Design and Control Laboratory
32
Natural convection in heterogeneous media
Two realizations of the porosity fields with 6
terms
Materials Process Design and Control Laboratory
33
Natural convection in heterogeneous media

• We have already obtained the KLE of the Gaussian
  random fields.

• Now We want to express the following
  non-polynomial functions in terms of polynomial
  chaos.

• We can evaluate the coefficients of the
  expansions for the three non-linear functions
  using LHS. Then we can sample the calculated GPCE
  expansion to obtain the PDF.
• The optimal order is determined such that the
  GPCE expansion can accurately represent the PDF
  of these non-linear functions.
• We compare the results with the Direct Sampling
  approach, where the PDF of the functions is
  obtained by sampling form the standard normal
  distribution, then calculate the realization of
  each sample and plot the PDF.

Materials Process Design and Control Laboratory
34
Natural convection in heterogeneous media
Probability
Probability
A second-order GPCE is enough to capture all of
the input uncertainties
Thus, we choose a 6-dimension second order GPCE
expansion to represent the solution process,
which give a total of 28 modes.
Probability
Materials Process Design and Control Laboratory
35
Natural convection in heterogeneous media
Mean
Materials Process Design and Control Laboratory
36
Natural convection in heterogeneous media
First order u velocity
Materials Process Design and Control Laboratory
37
Natural convection in heterogeneous media
First order v velocity
Materials Process Design and Control Laboratory
38
Natural convection in heterogeneous media
First order Temperature
Materials Process Design and Control Laboratory
39
Natural convection in heterogeneous media
Second order Modes
Materials Process Design and Control Laboratory
40
Natural convection in heterogeneous media
Probability
Probability
PDF at one point
velocity
velocity
Probability
Temperature
Materials Process Design and Control Laboratory
41
Natural convection in heterogeneous media
Standard deviation
Materials Process Design and Control Laboratory
42
Natural convection in heterogeneous media
Compare Mean
GPCE
Collocation
Monte Carlo
Materials Process Design and Control Laboratory
43
Natural convection in heterogeneous media
Compare standard deviation
GPCE
Collocation
Monte Carlo
Materials Process Design and Control Laboratory
44
3-D Natural convection in heterogeneous media
One realization of the random porosity random
field in 3-D
porosity slice of the xz plane at y0.5
KLE is truncated after two terms, so the
stochastic dimension is 2.
contour
The definition is the same as 2-D, except here we
consider a covariance function
porosity iso-surface
Materials Process Design and Control Laboratory
45
3-D Natural convection in heterogeneous media
Probability
Probability
Thus, we choose a 2-dimension third order GPCE
expansion to represent the solution process,
which give a total of 10 modes.
A second-order GPCE is not enough to capture all
of the input uncertainties. At least, we need a
third-order GPCE.
Probability
Materials Process Design and Control Laboratory
46
Mean slice of the xz plane at y0.5
Left GPCE
GPCE
Right Collocation
Collocation
Materials Process Design and Control Laboratory
47
Standard deviation slice of the xz plane at y0.5
Left GPCE
Right Collocation
GPCE
PDF at one point
Probability
Probability
Collocation
velocity
Temperature
Materials Process Design and Control Laboratory
48
Conclusions

• A second-order stabilized stochastic projection
  method was used to model uncertainty propagation
  in natural convection in random porous media.
• For the problems examined, the computation cost
  of the SSFEM and sparse grid simulation were
  similar. The MC simulations were significantly
  more expensive. In the SSFEM, it was shown that
  it is rather easy to identify dominant stochastic
  models in the solution and investigate how the
  uncertainty propagates from porosity to the
  velocity and temperature random fields.
• The key ingredient in the implementation of the
  algorithms presented here includes the
  development of a stochastic modeling library
  based on the GPCE formulation. This library
  includes computation tools for the implementation
  of the Askey polynomials, a parallel K-L
  expansion eigen solver, and a post-processing
  class for calculation of higher-order solution
  statistics such as standard deviation and
  probability density function.

Materials Process Design and Control Laboratory