An efficient surrogate model for emulation and physics extraction of large eddy simulations

1 Introduction

• High-fidelity simulations restrict practical design-space surveys; one injector simulation takes over 6 days on 200 CPU cores.The paper motivates emulation as a way to predict flows at new geometries using a small number of simulations.Source: §1 · p. 2
• The emulator uses flow-physics assumptions to reduce computation and improve prediction accuracy for massive simulation datasets.The paper also presents the assumptions as a way to make emulation practical for datasets spanning gigabytes or terabytes.Source: §1 · p. 2
• Existing flow-kriging methods face difficulty when spatial-grid dimensions vary across inputs, and the resulting flow data can be too large to process.For the injector application, the design range for length spans 20–100 mm.Source: §1 · p. 3
• Common POD extracts flow instabilities shared across design settings, addressing the difficulty of connecting structures across geometries.The method uses a physically justifiable linearity assumption about the spatial distribution of instability structures.Source: §1 · p. 3
• Domain-specific assumptions support uncertainty quantification, while a sparse co-kriging model learns dominant coupling mechanisms among flow variables.The resulting uncertainty can inform work on variable-disturbance propagation, and extracted physics can guide further investigation.Source: §1 · pp. 3–4

2.1 Injector design

• The injector design varies five geometric inputs—length, nozzle radius, inlet diameter, injection angle, and inlet-to-head-end distance—that affect film thickness and spreading angle.A larger injection angle is associated with thinner liquid films and smaller spreading angles.Source: §2.1 · p. 5

2 Injector schematic and large eddy simulations

• Large eddy simulation directly resolves large turbulent eddies and uses a model-based approach for small eddies.The section describes the simulations as the high-fidelity source of flow data for the emulator.Source: §2.2 · p. 6

2.2 Flow simulation

• The simulation dataset comprises 30 injector geometries, and each LES run takes 6 days on 200 CPU cores.The 30 runs are allocated across the design space using a maximum-projection design.Source: §2.2 · p. 6
• The 30 runs require nearly a hundred gigabytes of storage, while blind use of existing flow-kriging methods may require weeks for prediction.The paper presents physics-informed simplifications as a way to make prediction practical for this data scale.Source: §2.2 · pp. 6–7
• Accurate flow predictions can be achieved in hours by incorporating flow physics as simplifying assumptions, despite limited simulation runs.The paper states that 30 simulation runs can provide accurate predictions for this application.Source: §2.2 · pp. 6–7

3 Emulator model

• The emulator uses simplifying assumptions to support efficient, accurate flow prediction while extracting useful flow physics from data to inform future experiments.The authors state these objectives can be achieved despite limited simulation runs and complexities in flow data.Source: §3 · p. 7

3.1 Common POD

• CPOD extracts common flow structures across geometries by rescaling partitioned domains to a common grid under a linear-scaling assumption, justified for cold flows with similar Reynolds numbers.Its leading modes provide optimal data reduction and physically meaningful structures for emulation.Source: §3.1 · pp. 9–10

3.2 Model specification

• The model uses time-independent Gaussian processes with independent kriging at each time slice, avoiding temporal-correlation estimation when predictions between time steps are unnecessary and enabling parallel computation.The authors consider six flow variables and assume a separable covariance form for the coefficient vector.Source: §3.2 · pp. 11–12
• MA and independent model M0 have the same MMSE predictor when parameters are known, but MA improves uncertainty quantification and can predict derived flow functions differently.The paper demonstrates the difference for turbulent kinetic energy in Section 4.3.Source: §3.2 · p. 13

3.2.1 CPOD covariance matrix

• The covariance matrix T encodes physical couplings between flow modes, and selecting its correlations can rank significant couplings to guide further experiments.For example, injector pressure oscillations can induce velocity and density fluctuations.Source: §3.2.1 · p. 14
• With enough significant flow couplings, co-kriging model MA provides more accurate joint-prediction uncertainty quantification than independent model M0.The theorem indicates that M0 can have coverage below the desired 100(1 − α)% rate when it ignores these cross-correlations.Source: §3.2.1 · p. 16
• Imposing the element-wise L1 sparsity constraint ∥T−1∥1 ≤ γ, with small γ, makes nearly all entries of T−1 zero, permitting consistent estimation of the few significant correlations.The model permits different process variances for each coefficient, whereas the compared approach restricts all coefficients to equal variances.Source: §3.2.1 · pp. 16–17

3.3 Parameter estimation

• Maximum-likelihood estimation of μ, T, and τ is chosen over a Bayesian implementation for computational efficiency, so predictions can be generated faster than direct LES simulations requiring several days.The optimization uses a penalized negative log-likelihood.Source: §3.3 · p. 17
• The penalty λ controls sparsity in the selected correlations, and its tuning depends on whether predictive accuracy or correlation extraction is the goal.For predictive accuracy, the authors recommend cross-validation; for extracting correlations or using prior coupling information, they recommend setting λ to obtain a preset number of correlations.Source: §3.3 · pp. 17–18
• The blockwise coordinate descent algorithm alternates between optimizing the covariance matrix T and the Gaussian-process parameters µ and τ until convergence.The T update uses graphical LASSO; the µ and τ update uses a closed-form expression for µ and L-BFGS optimization for τ.Source: §3.3 · p. 18
• The BCD scheme converges to a stationary point, but it does not guarantee global optimization because the log-likelihood is non-convex in τ.The authors address this limitation by running parallel threads from different initial τ values and selecting the converged setting with the largest likelihood.Source: §3.3 · p. 19

4.1 Visualization and CPOD modes

• The leading two circumferential-velocity CPOD modes represent first- and second-harmonic longitudinal waves along the LOX film boundary, with an antinode in the second mode.The authors note that CPOD-mode interpretability enables extraction of physically meaningful couplings for further analysis.Source: Figure 6 · p. 20, §4.1 · p. 22
• Across n = 30 simulation cases, M = 10 and M = 45 modes captured 90% and 99% of total flow energy, respectively, versus around M = 20 for 99% in Zong and Yang (2008) for one geometry.The authors describe the CPOD as a reasonably compact representation and say it provides empirical evidence for the linearity assumption used for computation efficiency.Source: §4.1 · pp. 21–22

4.2 Emulation accuracy

• MRE stayed within 10% in the injector inlet, fluid transition region, and exit for nearly all fully developed-flow time steps (15–30 ms).The regions are evaluated because they are important for characterizing injector instability.Source: Figure 8 · p. 23, §4.2 · p. 25
• At the new geometry, predicted temperature flow closely matched the simulation in transition region, film thickness, spreading angle, surface waves, and downstream recirculation zones.The comparison demonstrates capture of key flow physics when known fluid-flow properties are incorporated as model assumptions.Source: Figure 7 · p. 23, §4.2 · p. 24
• Predicted and simulated probe spectra had nearly identical peaks at low and high frequencies, indicating that the emulator mimics important flow physics.The peaks can help identify physical properties of dominant instability structures.Source: §4.2 · p. 26

4.3 Uncertainty quantification

• Spatial uncertainty was lowest near the inlet and centerline and highest at three downstream gaseous cores, where recirculation vortices produce unstable flow.The pointwise confidence band covered realized prediction error and roughly matched its spatial distribution.Source: §4.3 · p. 27
• Temporal uncertainty at probe 1 was generally steady, with spikes around time-steps 300 and 800 linked respectively to central-vortex displacement and boundary development.The authors present uncertainty quantification as a way to extract flow physics without expensive simulations.Source: §4.3 · p. 27
• TKE prediction improved visually, and MA’s average coverage was 85.0% versus M0’s 73.8%, closer to the 90% nominal rate.The simulated TKE often dipped below M0’s lower confidence band.Source: §4.3 · pp. 28–29, Figure 13 · p. 29

4.4 Correlation extraction

• The extracted coupling graph links u1 and w1 through Bernoulli’s equation, w1 and P1 through centrifugal acceleration, and ρ1 and w2 through momentum conservation.The T1–ρ1 coupling reflects the equation of state and energy conservation.Source: §4.4 · pp. 29–30, Figure 14 · p. 30
• The CPOD matrix T serves both to improve uncertainty quantification and to extract flow couplings that can guide further experiments.The physical consistency of the selected couplings also reassures the authors about the accuracy of the estimated co-kriging matrix T.Source: §4.4 · p. 30

4.5 Computation time

• Despite a training dataset requiring nearly 100GB of storage, the emulator performs prediction, uncertainty quantification, and coupling extraction in slightly over an hour.After one-time CPOD extraction and parameter estimation, it can predict hundreds of new settings within a day.Source: §4.5 · p. 31
• The reported physical insights can guide further investigations.Source: §4.5 · p. 31

5 Conclusions and future work

• The method incorporates elicited flow properties as model assumptions and uses CPOD to extract common instabilities across the design space.The conclusion also describes time-independent emulators at each simulation time-step and a sparse covariance matrix T for significant flow-variable couplings.Source: §5 · p. 31
• The emulator accurately predicts flows and captures injector-performance metrics, using a fraction of the time required by simulation.The authors say this efficiency can support exploration of the desired design space and timely decisions.Source: §5 · pp. 31–32
• The model provides spatio-temporal uncertainty quantification and extracts significant couplings between flow instabilities.The authors describe the extracted couplings as useful for further analysis.Source: §5 · p. 32
• The authors propose pattern recognition within the GP kriging framework to identify instability structures that scale non-linearly across geometries and predict them at new settings.They identify computational efficiency as the key hurdle and mention treed GP and local GP models as options.Source: §5 · p. 32
• Future work includes evaluating whether MaxPro combined with minimax coverage can improve design performance.The authors refer to a recently proposed design in Mak and Joseph (2017).Source: §5 · p. 32
• The authors propose feeding emulated-flow uncertainty through an acoustics solver to evaluate new injector geometries, with representative points as a faster option.Each solver evaluation can be time-intensive; support points in Mak and Joseph (2016) are suggested for reducing the uncertainty propagation task.Source: §5 · p. 32

A Computing the CPOD expansion

• A common spatial domain is needed to extract common instabilities across injector geometries with different grid points.The appendix presents a physically justifiable method for obtaining the domain, then uses it to compute the CPOD expansion.Source: A Computing the CPOD expansion · p. 34

A.1 Common grid

• The procedure selects the simulation grid with the most points as the common reference grid.Each simulation grid is partitioned into four regions: injector head-end to inlet, inlet to nozzle exit, and the top and bottom downstream regions.Source: §A.1 · p. 34
• The four-part partition makes the linearity assumption physically justifiable.The partition separates injector and downstream regions before rescaling.Source: §A.1 · p. 34
• Each partition is linearly rescaled to the common grid using L, Rn, and ∆L, then flow data are interpolated onto its grid points.The implementation uses inverse-distance-weighting interpolation with 10 nearest neighbours.Source: §A.1 · p. 35

A.2 POD expansion

• The CPOD modes are computed from the common-grid snapshots and normalized to ensure orthonormality.The modes support extraction of common flow instabilities.Source: §A.2 · pp. 35–36
• The expansion indexes all combinations of n design settings and T time-steps as N = nT snapshots, each represented on common grid points.The inner-product matrix Q is defined over these snapshots.Source: §A.2 · p. 35
• The method computes eigenvectors of Q, whose eigenvalues rank the modes, but full eigendecomposition requires O(N^3) work.The implicitly restarted Arnoldi method is used to efficiently approximate leading eigenvalues and eigenvectors.Source: §A.2 · pp. 35–36
• Snapshot coefficients are derived for each design setting and time-step, then truncated at Kr < N modes to decompose the flow.The resulting decomposition is asserted in equation (3).Source: §A.2 · pp. 36–37

• Each BCD cycle applies S = A^M ◦ B^N, with M and N counting graphical LASSO iterations and line searches, respectively.Source: B Proof of Theorem 2 · p. 37
• The BCD iterates satisfy lim m→∞ Θm = Θ* ∈ Γ, where Γ contains points with ∇lλ(Θ)=0.Source: B Proof of Theorem 2 · p. 37
• The proof verifies three convergence conditions, including compact containment of the initial iterate and closedness of S outside Γ.Source: B Proof of Theorem 2 · pp. 37–38
• S preserves lλ on Γ and strictly decreases it off Γ because A and B decrease lλ for non-stationary states.Source: B Proof of Theorem 2 · p. 38

• Given simulation data D, the conditional flow vector y|D is Gaussian with mean prediction ŷ and covariance Φ(x,t).Source: C Proof of Theorem 3 · pp. 38–39
• After eigendecomposing Φ(t), κ(x,t) has a distribution that is a sum of weighted non-central chi-squared distributions.Source: C Proof of Theorem 3 · p. 39
• Pointwise confidence intervals for κ(x,t) use Liu et al. (2009)'s methods, implemented through the R package CompQuadForm.Source: C Proof of Theorem 3 · p. 40