Modeling diffusion in heterogeneous media: Data driven microstructure reconstruction models, stochastic collocation and

1. 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@cornell.edu http://mpdc.mae.cornell.edu * Work supported by AFOSR/Computational Mathematics

2. 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

3. 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

4. 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

5. 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. 1. PROPERTY EXTRACTION

7. 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

8. 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

9. 2. DATA DRIVEN MODELS FOR MICROSTRUCTURE RECONSTRUCTION

10. 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, 2π) Distributed according to where

11. 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

12. 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

13. 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

14. 3. REDUCED MODEL OF THE TOPOLOGICAL DESCRIPTOR

15. 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

16. REDUCED MODEL FOR THE STRUCTURE M microstructure images of nxnxn pixels each The microstructures are represented as vectors Ii i=1,..,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 +

17. 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

18. 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

19. 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

20. 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

21. 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

22. 4. SOLUTION TO THE STOCHASTIC PARTIAL DIFFERENTIAL EQUATION

23. SPDE Definition Governing equation for thermal diffusion Uncertainty comes in as the random material properties, which depend on the topology of the microstructure The (N+d) dimensional problem (N stochastic dimensions+ d spatial dimensions) is represented as The number of stochastic dimensions is usually large ~ 10-20

24. 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

25. 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

26. 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

27. 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

28. 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 Increasing the order of interpolation increases the number of points sampled

29. 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

30. 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

31. 5. SOLUTION TO THE DETERMINISTIC PARTIAL DIFFERENTIAL EQUATION

32. 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.

33. 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:

34. 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.

35. 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

36. 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

37. 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

38. 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

39. 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

40. 6. NUMERICAL EXAMPLE

41. MICROSTRUCTURE RECONSTRUCTION Experimental statistics Experimental image 3D microstructure GRF statistics

42. 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

43. 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 αAg/αW = 2.5 T= -0.5 T= 0.5 Left wall maintained at -0.5 Right wall maintained at +0.5 All other surfaces insulated

44. 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

45. FIRST ORDER STATISTICS: MEAN TEMPERATURE e c d b f a g

46. FIRST ORDER STATISTICS: HIGHER ORDER MOMENTS d c b e a f

47. SECOND ORDER STATISTICS: MEAN TEMPERATURE e c d b f g a

48. SECOND ORDER STATISTICS: HIGHER ORDER MOMENTS d c b e a f

49. 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 PDF’s at a point due to the application of first and second order constraints

50. 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