Bayesian Inference For The Calibration Of DSMC Parameters

1. Bayesian Inference For The Calibration Of DSMC Parameters James S. Strand and David B. Goldstein The University of Texas at Austin Sponsored by the Department of Energy through the PSAAP Program Computational Fluid Physics Laboratory Predictive Engineering and Computational Sciences

2. Motivation – DSMC Parameters • The DSMC model includes many parameters related to gas dynamics at the molecular level, such as: • Elastic collision cross-sections. • Vibrational and rotational excitation probabilities. • Reaction cross-sections. • Sticking coefficients and catalytic efficiencies for gas-surface interactions. • …etc.

3. DSMC Parameters • In many cases the precise values of some of these parameters are not known. • Parameter values often cannot be directly measured, instead they must be inferred from experimental results. • By necessity, parameters must often be used in regimes far from where their values were determined. • More precise values for important parameters would lead to better simulation of the physics, and thus to better predictive capability for the DSMC method. The ultimate goal of this work is to use experimental data to calibrate important DSMC parameters.

4. DSMC Method • Direct Simulation Monte Carlo (DSMC) is a particle based simulation method. • Simulated particles represent large numbers of real particles. • Particles move and undergo collisions with other particles. • Can be used in highly non-equilibrium flowfields (such as strong shock waves). • Can model thermochemistry on a more detailed level than most CFD codes.

5. DSMC Method Performed on First Time Step Initialize Create Move Performed on Every Time Step Index Collide Performed on Selected Time Steps Sample

6. DSMC Collisions DSMC collisions are performed statistically. Pairs of molecules are randomly selected from within a cell, a collision probability is calculated for each selected pair, and a random number draw determines whether or not the pair actually collides.

7. DSMC Collisions DSMC collisions are performed statistically. Pairs of molecules are randomly selected from within a cell, a collision probability is calculated for each selected pair, and a random number draw determines whether or not the pair actually collides.

8. DSMC Collisions DSMC collisions are performed statistically. Pairs of molecules are randomly selected from within a cell, a collision probability is calculated for each selected pair, and a random number draw determines whether or not the pair actually collides.

9. Selection Possibilities

10. DSMC Collisions If the pair does not collide, their properties are left unchanged and a new pair is selected. If the pair does collide, the properties of both colliding particles are immediately adjusted to their calculated post-collision values, and then a new pair is selected. The number of selections in a given cell is a function of the overall simulation conditions and of the conditions within that cell at that time step.

11. DSMC Collisions If the pair does not collide, their properties are left unchanged and a new pair is selected. If the pair does collide, the properties of both colliding particles are immediately adjusted to their calculated post-collision values, and then a new pair is selected. The number of selections in a given cell is a function of the overall simulation conditions and of the conditions within that cell at that time step.

12. Variable Hard Sphere Model The VHS model allows the collision cross-section to be dependent on relative speed. There are two relevant parameters for the VHS model, dref and ω. Both of these parameters are usually determined based on viscosity data.

13. Numerical Methods – DSMC Code • Our DSMC code can model flows with rotational and vibrational excitation and relaxation, as well as five-species air chemistry, including dissociation, exchange, and recombination reactions. • Larsen-Borgnakke model is used for redistribution between rotational, translational, and vibrational modes during inelastic collisions. • TCE model allows cross-sections for chemical reactions to be derived from Arrhenius parameters.

14. Internal Modes • Rotation is assumed to be fully excited. • Each particle has its own value of rotational energy, and this variable is continuously distributed. • Vibrational levels are quantized. • Each particle has its own vibrational level, which is associated with a certain vibrational energy based on the simple harmonic oscillator model. • Relevant parameters are ZR and ZV, the rotational and vibrational collision numbers. • ZR= 1/ΛR, where ΛR is the probability of the rotational energy of a given molecule being redistributed during a given collision. • ZV = 1/ΛV

15. Chemistry Implementation • Reaction cross-sections based on Arrhenius rates • TCE model allows determination of reaction cross-sections from Arrhenius parameters. σRand σTare the reaction and total cross-sections, respectively , the average number of internal degrees of freedom which contribute to the collision energy. is the temperature-viscosity exponent for VHS collisions between type A and type B particles k is the Boltzmann constant, mr is the reduced mass of particles A and B, Ec is the collision energy, and Γ() is the gamma function.

16. Reactions R. Gupta, J. Yos, and R. Thompson, NASA Technical Memorandum 101528, 1989.

17. Collision Rates

18. Reaction Rates – Dissociation

19. Reaction Rates – Recombination

20. Reaction Rates – Exchange

21. Parallelization • DSMC: • MPI parallel. • Ensemble averaging to reduce stochastic noise for 0-D relaxations. • Adaptive load balancing for 1-D shock simulations. • Sensitivity Analysis and MCMC: • MPI Parallel • Separate processor groups call DSMC subroutine independently.

22. Sensitivity Analysis - Overview • In the current context, the goal of sensitivity analysis is to determine which parameters most strongly affect a given quantity of interest (QoI). • Only parameters to which a given QoI is sensitive will be informed by calibrations based on data for that QoI.

23. Sensitivity Analysis Parameters Throughout the sensitivity analysis, the ratio of forward to backward rate for a given reaction is kept constant, since these ratios should be fixed by the equilibrium constant.

24. Scenario: 1-D Shock • Shock speed is ~8000 m/s, M∞≈ 23. • Upstream number density = 3.22×1021#/m3. • Upstream composition by volume: 79% N2, 21% O2. • Upstream temperature = 300 K.

25. 1D Shock Simulation • We require the ability to simulate a 1-D shock without knowing the post-shock conditions a priori. • To this end, we simulate an unsteady 1-D shock, and make use of a sampling region which moves with the shock.

26. Results: Nominal Parameter Values

27. Quantity of Interest (QoI) We cannot yet simulate NASA EASTor other shock tube results, so we must choose a temporary, surrogate QoI. ? J. Grinstead, M. Wilder, J. Olejniczak, D. Bogdanoff, G. Allen, and K. Danf, AIAA Paper 2008-1244, 2008.

28. Sensitivity Analysis - Scalar vs. Vector QoI The density of NO (our QoI for this work) is a vector, with values over a range of points in space.

29. Sensitivity Analysis - Methods • Two measures for sensitivity were used in this work. • Pearson correlation coefficients: • The mutual information: • Both measures involve global sensitivity analysis based on a Monte Carlo sampling of the parameter space, and thus the same datasets can be • used to obtain both measures.

30. Sensitivity Analysis: Correlation Coefficient

31. Sensitivity Analysis: Mutual Information

32. Sensitivity Analysis: Mutual Information

33. Sensitivity Analysis: Mutual Information

34. Sensitivity Analysis: Mutual Information Actual joint PDF of θ1and the QoI, from a Monte Carlo sampling of the parameter space. Actual joint PDF of θ1and the QoI, from a Monte Carlo sampling of the parameter space. Kullback-Leibler divergence QoI Hypothetical joint PDF for case where the QoI is indepenent of θ1. Hypothetical joint PDF for case where the QoI is indepenent of θ1. θ1

35. Sensitivities vs. X

36. r2 vs. Mutual Information

37. Sensitivities vs. X for ρNOas QoI

38. Variance of ρNOvs. X

39. Variance Weighted Sensitivities

40. Overall Sensitivities

41. Synthetic Data Calibrations: 0-D Relaxation • Due to computational constraints, a 0-D relaxation from an initial high-temperature state was used for the synthetic data calibrations performed thus far. • 0-D box is initialized with 79% N2, 21% O2. • Initial bulk number density = 1.0×1023#/m3. • Initial bulk translational temperature = ~50,000 K. • Initial bulk rotational and vibrational temperatures are both 300 K. • Scenario is a 0-D substitute for a hypersonic shock at ~8 km/s. • Assumption that the translational modes equilibrate much faster than the internal modes. • Sensitivity analysis of the type discussed earlier was used to identify parameters for calibration.

42. Synthetic Data • We once again use ρNO as our QoI.

43. MCMC Method - Overview • Markov Chain Monte Carlo (MCMC) is a method which solves the statistical inverse problem in order to calibrate parameters with respect to a set of data. • The likelihood of a given set of parameters is calculated based on the mismatch between the data and the simulation results for that set of parameters. • One or more chains explore the parameter space, moving towards regions of higher likelihood. • Candidate positions are drawn from a multi-dimensional Gaussian proposal distribution centered at the current chain position. The covariance matrix of this Gaussian controls the average distance (in parameter space) that the chain moves in one step.

44. MCMC Method - Likelihood

45. MCMC Method – Metropolis-Hastings Establish boundaries for parameter space Candidate Accepted Candidate Rejected Accept or reject candidate position based on a random number draw Select initial position Current position remains unchanged. Candidate position becomes current position Run simulation at current position Likelihoodcandidate < Likelihoodcurrent Run simulation for candidate position parameters, and calculate likelihood Calculate likelihood for the current position Draw new candidate position Likelihoodcandidate > Likelihoodcurrent Candidate position is accepted, and becomes the current chain position Candidate automatically accepted

46. MCMC Method - Improvements • We use the PECOS-developed code QUESO, which implements two major additions to the basic Metropolis-Hastings algorithm. Both of these features can help improve the convergence of the method. • Delayed Rejection: • When an initial candidate position is rejected, a second candidate position is generated based on a scaled proposal covariance matrix. • Adaptive Metropolis: • At periodic intervals, the proposal covariance matrix is updated based on the calculated covariance of the previously accepted chain positions. H. Haario, M. Laine, A. Mira, E. Saksman, Statistics and Computing, 16, 339-354 (2006).

47. Chain Progression

48. Post-Calibration PDFs

49. Conclusions • Global, Monte Carlo based sensitivity analysis can provide a great deal of insight into how various parameters affect a given QoI. • A great deal of computer power is required to perform this type of statistical analysis for DSMC shock simulations. • 5,000 shocks were run for this sensitivity analysis. • Required a total of ~320,000 CPU hours. • MCMC can be used to calibrate at least some DSMC parameters based on synthetic data for a 0-D relaxation. • MCMC allows for propagation of uncertainty from the data to the final parameter PDFs.

50. Future Work • Synthetic data calibrations for a 1-D shock with the current code. • Upgrade the code to allow modeling of ionization and electronic excitation. • Couple the code with a radiation solver. • Sensitivity analysis for a 1-D shock with the additional physics included. • Synthetic data calibrations with the upgraded code. • Calibrations with real data from EAST or similar facility.