1 / 50

Modeling diffusion in heterogeneous media Data

driven microstructure reconstruction models,

stochastic collocation and the variational

multiscale method Nicholas Zabaras and Baskar

Ganapathysubramanian Materials Process Design and

Control Laboratory Sibley School of Mechanical

and Aerospace Engineering Cornell

University Ithaca, NY 14853-3801 zabaras_at_cornell.

edu http//mpdc.mae.cornell.edu Work

supported by AFOSR/Computational Mathematics

TRANSPORT IN HETEROGENEOUS MEDIA

- Thermal and fluid transport in heterogeneous

media are ubiquitous - 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 error

propagate through them?

- To apply physical processes on these

heterogeneous systems - worst case scenarios
- variations on physical properties

ISSUES WITH INVESTIGATION TRANSPORT IN

HETEROGENEOUS MEDIA

- Some critical issues have to be resolved to

achieve realistic results. - Multiple length scale variations in the material

properties of the heterogeneous medium - The essentially statistical nature of information

available about the media - Presence of uncertainty in the system and

properties

Only some statistical features can be extracted

PROBLEM OF INTEREST

Interested in modeling diffusion through

heterogeneous random media

Aim To develop procedure to predict statistics

of properties of heterogeneous materials

undergoing diffusion based transport

- Should account for the multi length scale

variations in thermal properties - Account for the uncertainties in the topology of

the heterogeneous media

- What is given
- Realistically speaking, one usually has access to

a few experimental 2D images of the

microstructure. Statistics of the heterogeneous

microstructure can then be extracted from the

same. - This is our starting point

OVERVIEW OF METHODOLOGY

2. Microstructure reconstruction

1. Property extraction

Extract properties P1, P2, .. Pn, that the

structure satisfies. These properties are usually

statistical Volume fraction, 2 Poit correlation,

auto correlation

Reconstruct realizations of the structure

satisfying the properties. Monte Carlo, Gaussian

Random Fields, Stochastic optimization ect

Construct a reduced stochastic model from the

data. This model must be able to approximate the

class of structures. KL expansions, FFT and other

transforms, Autoregressive models, ARMA models

Solve the heterogeneous property problem in the

reduced stochastic space for computing property

variations. Collocation schemes VMS

4. Stochastic collocation Variational

multiscale method

3. Reduced model

1. PROPERTY EXTRACTION

IMAGE PROCESSING

Reconstruction of well characterized

material Tungsten-Silver composite1 Produced by

infiltrating porous tungsten solid with molten

silver

640x640 pixels 198 Âµm x 198 Âµm

1. S. Umekawa, R. Kotfila, O. D. Sherby, Elastic

properties of a tungsten-silver composite above

and below the melting point of silver, J. Mech.

Phys. Solids 13 (1965) 229-230

PROPERTY EXTRACTION

First order statistics Volume fraction 0.2

Second order statistics 2 pt correlation

Digitized two phase microstructure image White

phase- W Black phase- Ag Simple matrix operations

to extract image statistics

2. DATA DRIVEN MODELS FOR MICROSTRUCTURE

RECONSTRUCTION

MICROSTRUCTURE RECONSTRUCTION

Statistical information available- First and

second order statistics Reconstruct Three

dimensional microstructures that satisfy these

experimental statistical relations

GAUSSIAN RANDOM FIELDS GRF- model interfaces as

level cuts of a function Build a function y(r).

Model microstructure is given by level cuts of

this function. y(r) has a field-field correlation

given by g(r) If this function is known, y(r) can

be constructed as

Uniformly distributed over the unit sphere

Uniformly distributed over 0, 2p)

Distributed according to where

MICROSTRUCTURE RECONSTRUCTION

- Relate experimental properties to
- Two phase microstructure, impose level cuts on

y(r). Phase 1 if - Relate to statistics
- first order statistics

where

second order statistics

Set , and For the Gaussian Random Field

to match experimental statistics

MICROSTRUCTURE RECONSTRUCTION FITTING THE GRF

PARAMETERS

Assume a simplified form for the far field

correlation function

Three parameters, ÃŸ is the correlation length, d

is the domain length and rc is the cutoff length

Use least square minimization to find optimal fit

3D MICROSTRUCTURE RECONSTRUCTION

20 Âµm x 20 Âµm x 20 Âµm

64x64x64 pixel

40 Âµm x 40 Âµm x 40 Âµm

128x128x128 pixel

200 Âµm x 200 Âµm

3. REDUCED MODEL OF THE TOPOLOGICAL DESCRIPTOR

WHY A REDUCED MODEL?

The reconstruction procedure gives a large set of

3D microstructures The topology of the

reconstructured microstructures are all

different All these structures satisfy the

experimental statistical relations These

microstructures belong to a very large (possibly)

infinite dimensional space. These topological

variations are the inputs to the stochastic

problem The necessity of model reduction arises

Model reduction techniques Most commonly used

technique in this context is Principle Component

Analysis Compute the eigen values of the dataset

of microstructures

REDUCED MODEL FOR THE STRUCTURE

M microstructure images of nxnxn pixels each The

microstructures are represented as vectors Ii

i1,..,M The eigenvectors of the n3xn3 covariance

matrix are computed The first N eigenimages are

chosen to represent the microstructures

Represent any microstructure as a linear

combination of the eigenimages

I Iavg I1a1 I2a2 I3a3 Inan

an

..

a2

a1

REDUCED MODEL FOR THE STRUCTURE CONSTRAINTS

Let I be an arbitrary microstructure satisfying

the experimental statistical correlations The PCA

method provides a unique representation of the

image That is, the PCA provides a function

The function is injective but nor

surjective Every image has a unique

mapping But every point

need not define an image in

Construct the subspace of allowable n-tuples

CONSTRUCTING THE REDUCED SUBSPACE H

Image I belongs to the class of structures? It

must satisfy certain conditions a) Its volume

fraction must equal the specified volume

fraction b) Volume fraction at every pixel must

be between 0 and 1 c) It should satisfy the

given two point correlation Thus the n tuple

(a1,a2,..,an) must further satisfy some

constraints. Enforce these constraints

sequentially

1. Pixel based constraints

Microstructures represented as discrete images.

Pixels have bounds This results in 2n3 inequality

constraints

CONSTRUCTING THE REDUCED SUBSPACE H

2. First order constraints

The Microstructure must satisfy the experimental

volume fraction

This results in one linear equality constraint on

the n-tuple

3. Second order constraints

The Microstructure must satisfy the experimental

two point correlation. This results in a set of

quadratic equality constraints

This can be written as

SEQUENTIAL CONSTRUCTION OF THE SUBSPACE

Computational complexity Pixel based constraints

first order constraints result in a simple

convex hull problem Enforcing second order

constraints becomes a problem in quadratic

programming Sequential construction of the

subspace First enforce first order statistics, On

this reduced subspace, enforce second order

statistics Example for a three dimensional space

3 eigen images

THE REDUCED MODEL

The sequential contraction procedure a subspace

H, such that all n-tuples from this space result

acceptable microstructures

H represents the space of coefficients that map

to allowable microstructures. Since H is a plane

in N dimensional space, we call this the

material plane

Since each of the microstructures in the

material plane satisfies all required

statistical properties, they are equally

probable. This observation provides a way to

construct the stochastic model for the allowable

microstructures

Define such that

This is our reduced stochastic model of the

random topology of the microstructure class

4. SOLUTION TO THE STOCHASTIC PARTIAL

DIFFERENTIAL EQUATION

SPDE Definition

Governing equation for thermal diffusion

Uncertainty comes in as the random material

properties, which depend on the topology of the

microstructure

The (Nd) dimensional problem (N stochastic

dimensions d spatial dimensions) is represented

as

The number of stochastic dimensions is usually

large 10-20

UNCERTAINTY ANALYSIS TECHNIQUES

- Monte-Carlo Simple to implement,

computationally expensive - Perturbation, Neumann expansions Limited to

small fluctuations, tedious for higher order

statistics

- Spectral stochastic uncertainty representation

Basis in probability and functional analysis, Can

address second order stochastic processes, Can

handle large fluctuations, derivations are

general - Stochastic collocation Results in decoupled

equations

COLLOCATION TECHNIQUES

Spectral Galerkin method Spatial domain is

approximated using a finite element

discretization Stochastic domain is

approximated using a spectral element

discretization

Coupled equations

Decoupled equations

Collocation method Spatial domain is

approximated using a finite element

discretization Stochastic domain is

approximated using

multidimensional interpolating functions

DECOUPLED EQUATIONS IN STOCHASTIC SPACE

Simple interpolation Consider the function We

evaluate it at a set of points The approximate

interpolated polynomial representation for the

function is Where Here, Lk are the Lagrange

polynomials Once the interpolation function has

been constructed, the function value at any point

yi is just Considering the given natural

diffusion system One can construct the

stochastic solution by solving at the M

deterministic points

SMOLYAK ALGORITHM

LET OUR BASIC 1D INTERPOLATION SCHEME BE

SUMMARIZED AS

IN MULTIPLE DIMENSIONS, THIS CAN BE WRITTEN AS

TO REDUCE THE NUMBER OF SUPPORT NODES WHILE

MAINTAINING ACCURACY WITHIN A LOGARITHMIC FACTOR,

WE USE SMOLYAK METHOD

IDEA IS TO CONSTRUCT AN EXPANDING SUBSPACE OF

COLLOCATION POINTS THAT CAN REPRESENT

PROGRESSIVELY HIGHER ORDER POLYNOMIALS IN

MULTIPLE DIMENSIONS A FEW FAMOUS SPARSE

QUADRATURE SCHEMES ARE AS FOLLOWS CLENSHAW

CURTIS SCHEME, MAXIMUM-NORM BASED SPARSE GRID AND

CHEBYSHEV-GAUSS SCHEME

SMOLYAK ALGORITHM

Extensively used in statistical

mechanics Provides a way to construct

interpolation functions based on minimal number

of points Univariate interpolations to

multivariate interpolations

Uni-variate interpolation

Multi-variate interpolation

Smolyak interpolation

Accuracy the same as tensor product Within

logarithmic constant

D 10

ORDER SC FE

3 1581 8000

4 8801 40000

5 41625 100000

Increasing the order of interpolation increases

the number of points sampled

SMOLYAK ALGORITHM REDUCTION IN POINTS

For 2D interpolation using Chebyshev nodes Left

Full tensor product interpolation uses 256

points Right Sparse grid collocation used 45

points to generate interpolant with comparable

accuracy

D 10

Results in multiple orders of magnitude

reduction in the number of points to sample

ORDER SC FE

3 1581 1000

4 8801 10000

5 41625 100000

SPARSE GRID COLLOCATION METHOD implementation

Solution Methodology

PREPROCESSING Compute list of collocation points

based on number of stochastic dimensions, N and

level of interpolation, q Compute the weighted

integrals of all the interpolations functions

across the stochastic space (wi)

Use any validated deterministic solution

procedure. Completely non intrusive

Solve the deterministic problem defined by each

set of collocated points

POSTPROCESSING Compute moments and other

statistics with simple operations of the

deterministic data at the collocated points and

the preprocessed list of weights

Std deviation of temperature Natural convection

5. SOLUTION TO THE DETERMINISTIC PARTIAL

DIFFERENTIAL EQUATION

THE NECESSITY FOR VARIATIONAL MULTISCALE METHODS

The collocation method reduces the stochastic

problem to the solution of a set of deterministic

equations

- These deterministic problems correspond to

solving the thermal diffusion problem on a set of

unique microstructures. - These heterogeneous

microstructure realizations exhibit property

variations at a much smaller scale compared to

the size of the computational domain - Performing

a fully-resolved calculation on these

microstructures becomes computationally

expensive. - Consider a computational scheme that

involves solving for a coarse-solution while

capturing the effects of the fine scale on the

coarse solution.

ADDITIVE SCALE DECOMPOSITION

The variation form of the diffusion equation can

be written as

- Assume that the solution can be decomposed into

two scales - Coarse resolvable scale
- Fine irresolvable (but modeled) scale

The variation form of the diffusion equation

decomposes into

SUB GRID MODELLING

Further decompose fine scale solution into two

parts

Particular solution

Homogeneous solution

- The solution component incorporates the

entire coarse scale solution information and

has no dependence on the coarse scale solution.

- The dynamics of is driven by the

projection of the source term onto the

subgrid scale function space.

SUB GRID MODELLING II

piecewise polynomial finite element

representation for the coarse solution inside a

coarse element Similar representation for the

fine scale. Move problem from computing values at

finest resolution to computing the shape function

at the finest resolution

Substitute into fine scale variational equation

SUB GRID MODELLING III

Without loss of generality, we can assume the

following representation for the coarse scale

nodal solutions

A very general representation that incorporates

several well known time integration schemes

Substituting this form for the coarse and fine

scale solutions into the fine scale variational

forms gives

This is valid for all possible combinations of u.

It follows that each of the quantities in the

brackets above must equal 0

SUB GRID MODELLING IV

This gives the variational form for the sub-grid

basis functions

The strong form for the fine-scale basis function

is then given by

The solution of the fine scale evolution equation

can then be input into the coarse scale solution

to get the coarse scale evolution equation

VERIFICATION OF THE VMS FORMULATION

Reconstructed VMS solutions

Coarse scale VMS solutions

(a) Fully resolved FEM solution

(d)

(e)

(c)

(b)

Increasing coarse element size

OVERVIEW OF METHODOLOGY

2. Microstructure reconstruction

1. Property extraction

Extract properties P1, P2, .. Pn, that the

structure satisfies. These properties are usually

statistical Volume fraction, 2 Poit correlation,

auto correlation

Reconstruct realizations of the structure

satisfying the properties. Monte Carlo, Gaussian

Random Fields, Stochastic optimization ect

Construct a reduced stochastic model from the

data. This model must be able to approximate the

class of structures. KL expansions, FFT and other

transforms, Autoregressive models, ARMA models

Solve the heterogeneous property problem in the

reduced stochastic space for computing property

variations. Collocation schemes VMS

4. Stochastic collocation Variational

multiscale method

3. Reduced model

6. NUMERICAL EXAMPLE

MICROSTRUCTURE RECONSTRUCTION

Experimental statistics

Experimental image

3D microstructure

GRF statistics

MODEL REDUCTION

Principal component analysis

Constructing the reduced subspace and the

stochastic model

- Enforcing the pixel based bounds and the linear

equality constraint (of volume fraction) was

developed as a convex hull problem. A primal-dual

polytope method was employed to construct the set

of vertices. - Enforcing the second order constraints was

performed through the quadratic programming tools

in the optimization toolbox in Matlab. - Two separate cases are considered in this
- example. In the first case, only the first-order

constraints (volume fraction) are used to

reconstruct the subspace H. In the second case,

both first-order as well as second-order

constraints (volume fraction and two-point

correlation) are used to construct the subspace H.

First 9 eigen values from the spectrum chosen

PHYSICAL PROBLEM UNDER CONSIDERATION

Structure size 40x40x40 Âµm Tungsten Silver

Matrix Heterogeneous property is the thermal

diffusivity. Tungsten ? 19250 kg/m3 k 174

W/mK c 130 J/kgK Silver ? 10490 kg/m3 k

430 W/mK c 235 J/kgK Diffusivity ratio aAg/aW

2.5

T -0.5

T 0.5

Left wall maintained at -0.5 Right wall

maintained at 0.5 All other surfaces insulated

COMPUTATIONAL DETAILS

The construction of the stochastic solution

through sparse grid collocation level 5

interpolation scheme used Number of deterministic

problems solved 15713

Computational domain of each deterministic

problem 128x128x128 pixels

Each deterministic problem solution solved on a

8 8 8 coarse element grid (uniform hexahedral

elements) with each coarse element having 16 16

16 fine-scale elements.

The solution of each deterministic VMS problem

about 34 minutes, In comparison, a fully-resolved

fine scale FEM solution took nearly 40 hours.

Computational platform 40 nodes on local Linux

cluster Total time 56 hours

FIRST ORDER STATISTICS MEAN TEMPERATURE

e

c

d

b

f

a

g

FIRST ORDER STATISTICS HIGHER ORDER MOMENTS

d

c

b

e

a

f

SECOND ORDER STATISTICS MEAN TEMPERATURE

e

c

d

b

f

g

a

SECOND ORDER STATISTICS HIGHER ORDER MOMENTS

d

c

b

e

a

f

CONCLUSIONS

A new model for modeling diffusion in random

two-phase media. A general methodology was

presented for constructing a reduced-order

microstructure model for use as random input in

the solution of stochastic partial differential

equations governing physical processes The twin

problems of uncertainty and multi length scale

variations are decoupled and comprehensively

solved Scope of further research Using more

sophisticated model reduction techniques to build

the reduced-order microstructure model, Extending

the methodology to arbitrary types of

microstructures as well as developing models of

advection-diffusion in random heterogeneous

media.

Comparison of temperature PDFs at a point due to

the application of first and second order

constraints

REFERENCES

- B. Ganapathysubramanian and N. Zabaras, "Sparse

grid collocation methods for stochastic natural

convection problems", Journal of Computational

Physics, in press - B. Ganapathysubramanian and N. Zabaras,

"Modelling diffusion in random heterogeneous

media Data-driven models, stochastic collocation

and the variational multi-scale method", Journal

of Computational Physics, submitted - S. Sankaran and N. Zabaras, "Computing property

variability of polycrystals induced by grain size

and orientation uncertainties", Acta Materialia,

in press - B. Velamur Asokan and N. Zabaras, "A stochastic

variational multiscale method for diffusion in

heterogeneous random media", Journal of

Computational Physics, Vol. 218, pp. 654-676,

2006 - B. Velamur Asokan and N. Zabaras, "Using

stochastic analysis to capture unstable

equilibrium in natural convection", Journal of

Computational Physics, Vol. 208/1, pp. 134-153,

2005