1 Introduction↩︎

Non-intrusive reduced-order emulators substitute for repeated executions of expensive finite-element solvers by training cheap approximations on solver snapshots. This approach enables design studies and many-query applications that remain computationally intractable with the solver alone. When targeted at certification, the utility of such emulators depends on establishing rigorous error bounds. These bounds require evaluation against the correct reference. Rather than physical truth, the reference is the underlying numerical solver. This solver is itself an approximation of physical reality and possesses its own discretization error. Reducing the emulator error below this discretization threshold merely refines agreement with the numerical model beyond its own resolution, which lacks physical meaning.

Intrusive reduced-basis methods benefit from established a posteriori error estimators that access the governing operators to bound the reduced solution rigorously [1], [2]. In contrast, the non-intrusive setting operates strictly on input-output data. Because residual-based estimators are unavailable in this context, establishing a controlled, certification-grade error assessment against the solver’s numerical uncertainty remains an open challenge. Standard non-intrusive reduced-order models construct maps via regression of reduced coordinates or via operator inference [3]–[5]. Separate uncertainty and calibration frameworks map the discrepancy between models and physical systems [6] within the standard lexicon of verification and validation [7].

This study evaluates a high-consequence benchmark from aircraft braking. A rejected take-off converts the entire kinetic energy of the aircraft into thermal energy, causing the carbon-carbon brake disc to experience a peak temperature followed closely by peak thermoelastic stress [8]. These peak values dictate landing gear safety and brake reuse. Predicting them across the full operating envelope represents a primary design and certification requirement. High-fidelity finite-element models of the disc provide accurate simulations. Their computational cost makes dense parametric exploration impractical for multi-query studies [8]. This limitation motivates data-driven surrogates for virtual sensing of quantities that are hard to measure directly [9], prediction of the brake stress state from simulation data [10], and real-time structural assessment [11].

This paper presents a verification methodology for non-intrusive reduced-order emulators derived from high-fidelity solvers, demonstrated on the aircraft brake disc. The underlying emulator uses a standard pipeline coupling proper orthogonal decomposition (POD) and Gaussian process regression (GPR). The novelty lies in the verification framework. We introduce an a posteriori peak error budget that decomposes the signed emulator error, via an algebraically exact identity, into a reduced-basis truncation term and a regression term. The observed extremum error is then evaluated against a discretization floor derived from a mesh-convergence study of the reference model. Predicting peak values is a difficult task. The extremum is a non-smooth functional that amplifies local spatial errors, making peak verification more stringent than global field-wise norms. Our main finding shows that, on a held-out test set, the emulator reproduces the high-fidelity peak predictions within the solver’s own discretization uncertainty. We define this level of agreement as surrogate-to-code fidelity, which extends code verification concepts to surrogate models rather than asserting physical validation. Two secondary contributions support this framework. We perform concurrent non-intrusive predictions of thermal and mechanical peaks within a single workflow. We also introduce a probabilistic localization of spatial extrema zones, propagating Gaussian process posterior uncertainties directly onto the computational mesh, reported under an explicit appraisal of its added value.

The manuscript is organized as follows. Section 2 outlines the mathematical formulation and the peak quantities of interest. Section 3 defines the high-fidelity model, the dataset, and the discretization floor. Section 4 details the non-intrusive emulator pipeline. Section 5 derives the peak error budget and the level-one uncertainty metrics. Section 6 presents the results obtained on the held-out test set. Section 7 provides an ablation study, evaluates the scope of the linear reduction, and discusses model interpretability. Section 8 addresses limitations and positioning, and Section 9 concludes the paper.

2 Problem statement and quantities of interest at peak↩︎

Brake certification protocols rely on the peak thermal state of a brake disc to evaluate structural integrity. This maximum thermal state also guides operational decisions regarding brake reuse and aircraft turnaround times after an aborted take-off. A rejected take-off represents the most severe thermal loading case for sizing braking systems. Modern high-energy aircraft designs utilize carbon-carbon composite discs to exploit their high thermal capacity and minimal mass [8]. In these scenarios, the brakes convert the kinetic energy of the aircraft entirely into thermal energy, and this heat generation drives the subsequent mechanical response of the disc. This study focuses on this peak operating state. We develop a non-intrusive emulator to predict both the temperature field and the von Mises stress field of a carbon-carbon brake disc at their respective peaks using only the operating parameters as inputs. The model also provides predictive uncertainty estimates for these peak quantities relative to the high-fidelity solver predictions.

2.1 Parametric input space↩︎

We define each operating condition by a vector \(\mathbf{x} \in \mathbb{R}^{5}\) containing 5 scalar parameters: \[\mathbf{x} = \bigl(P,\; M,\; V_0,\; T_0,\; T_{\mathrm{amb}}\bigr), \label{eq:input-vector}\tag{1}\] where \(P\) represents the braking load in kilonewtons, \(M\) is the aircraft mass in tonnes, \(V_0\) is the initial speed during the rejected take-off in meters per second, \(T_0\) denotes the initial disc temperature in degrees Celsius, and \(T_{\mathrm{amb}}\) is the ambient air temperature in degrees Celsius. Section 3 details the admissible parameter ranges and the sampling strategy used to build the dataset.

2.2 Prediction targets↩︎

For each operating point, the high-fidelity model described in Section 3 computes the transient thermal and mechanical response over a fixed mesh containing \(N_{\mathrm{n}} = 60914{}\) nodes. This shared discretization guarantees direct nodal correspondence across all dataset samples. We denote the transient temperature and von Mises stress fields as \(\mathcal{T}(\mathbf{x},t) \in \mathbb{R}^{N_{\mathrm{n}}}\) and \(\mathcal{S}(\mathbf{x},t) \in \mathbb{R}^{N_{\mathrm{n}}}\). Rather than emulating the complete transient history, we extract each field at the specific time step where its global spatial maximum is reached: \[t_T(\mathbf{x}) = \arg\max_{t}\;\max_{1 \le i \le N_{\mathrm{n}}} \mathcal{T}_i(\mathbf{x},t), \qquad t_\sigma(\mathbf{x}) = \arg\max_{t}\;\max_{1 \le i \le N_{\mathrm{n}}} \mathcal{S}_i(\mathbf{x},t). \label{eq:peak-instants}\tag{2}\] Focusing on these discrete instances aligns with certification requirements, which prioritize peak thermal and mechanical loads over full transient histories. The resulting prediction targets are the spatial fields: \[T(\mathbf{x}) = \mathcal{T}\bigl(\mathbf{x}, t_T(\mathbf{x})\bigr) \in \mathbb{R}^{N_{\mathrm{n}}}, \qquad \sigma(\mathbf{x}) = \mathcal{S}\bigl(\mathbf{x}, t_\sigma(\mathbf{x})\bigr) \in \mathbb{R}^{N_{\mathrm{n}}}, \label{eq:target-fields}\tag{3}\] representing the nodal temperature (in degrees Celsius) and the von Mises stress (in megapascals) at their respective peak times. In our dataset, the mechanical peak consistently occurs at or after the thermal peak, yielding \(t_\sigma(\mathbf{x}) \ge t_T(\mathbf{x})\). This is an empirical property of the present dataset and is not imposed by the formulation. The emulator learns these two static states directly, bypassing trajectory modeling.

2.3 Certified quantities of interest↩︎

The target certified values are scalar output functionals derived from the predicted spatial fields, following established principles for parametrized problems [12]. For any nodal field \(u \in \mathbb{R}^{N_{\mathrm{n}}}\), we define the peak operator as: \[Q(u) = \max_{1 \le i \le N_{\mathrm{n}}} u_i , \label{eq:peak-functional}\tag{4}\] yielding the two primary quantities of interest, \(Q\bigl(T(\mathbf{x})\bigr)\) and \(Q\bigl(\sigma(\mathbf{x})\bigr)\). This functional \(Q\) represents a non-differentiable, nonlinear \(L^{\infty}\) extremum evaluated over the mesh nodes. Operating on discrete vectors makes \(Q\) a mesh-dependent approximation of the continuous spatial supremum, supporting our focus on a surrogate-to-code verification strategy. Standard reduced-basis output bounding focuses on linear functionals within coercive systems [12]. Since \(Q\) is a nonlinear maximum operator, we certify the peak values indirectly using a posteriori bounds on the spatial field error instead of direct analytical output bounds. The high-fidelity finite-element simulations in Section 3 serve as the reference. Verifying the emulator directly against solver outputs establishes a surrogate-to-code rather than a surrogate-to-physics certification. Section 5 develops the exact a posteriori peak error budget supporting this verification.

2.4 Spatial structure of the extremum↩︎

While the peak scalar \(Q(u)\) remains mathematically well-posed, the specific node where the maximum occurs lacks robustness as a physical descriptor. The high-fidelity fields in Figure 1 show that the extrema occupy a broad, nearly flat annular band on the disc friction surface. The maximum value varies smoothly across the operating space, whereas the argmax coordinates shift abruptly within this flat crest, preventing reliable node-level identification. This flat profile near the maximum justifies identifying the extremum location as a spatial region rather than a single nodal point. For any field \(u\) and tolerance parameter \(\alpha \in (0,1)\), we define this region as the super-level set: \[Z_\alpha(u) = \bigl\{\, i \in \{1,\dots,N_{\mathrm{n}}\} \;\mid\; u_i \ge (1-\alpha)\,Q(u) \,\bigr\}. \label{eq:zone}\tag{5}\] We develop the probabilistic prediction of this spatial zone as the secondary contribution of this study, with an explicit caveat on its added value in Section 7. The formal a posteriori peak error budget derived in Section 5 constitutes the primary methodology of our verification framework.

3 High-fidelity model, dataset, and the discretization floor↩︎

This section details the reference finite-element model used to train and test the emulator. We outline the dataset generation, the specific modeling assumptions bounding the surrogate-to-code claim, and the discretization floor that defines the limit of meaningful emulator accuracy.

3.1 Governing thermomechanical model↩︎

A rejected take-off converts the kinetic energy of the aircraft into heat at the disc friction interface, generating a transient temperature field that subsequently drives thermoelastic stress. This coupling is unidirectional, meaning temperature changes alter the stress state while mechanical work does not affect the thermal field. The disc occupies a three-dimensional domain \(\Omega \subset \mathbb{R}^{3}\).

The thermal field satisfies the transient heat equation. With density \(\rho\), specific heat \(c_p\), and orthotropic conductivity tensor \(\mathbf{k}\), the energy balance and Fourier’s law yield \[\rho\, c_p\, \frac{\partial T}{\partial t} = \nabla \cdot \bigl( \mathbf{k}\, \nabla T \bigr). \label{eq:heat}\tag{6}\] The carbon-carbon composite is orthotropic, with a diagonal conductivity tensor in the material axes, \(\mathbf{k} = \operatorname{diag}(k_1, k_2, k_3)\), reflecting differing in-plane and through-thickness transport. Frictional heating acts strictly as a boundary flux, meaning no volumetric source term enters 6 and no mechanical dissipation is added. We assume temperature-independent material properties over the analyzed operating range. The initial condition imposes a uniform temperature, \(T = T_0\) at \(t = 0\). Boundary conditions combine the frictional heat flux on contact surfaces with convective cooling on free surfaces, \[\begin{align} -\mathbf{k}\, \nabla T \cdot \mathbf{n} &= \dot{q}_0(t) && \text{(Neumann, contact faces)}, \\ -\mathbf{k}\, \nabla T \cdot \mathbf{n} &= h\,(T - T_{\mathrm{amb}}) && \text{(Robin, free faces)}. \end{align} \label{eq:thermal-bc}\tag{7}\] The transient flux magnitude \(\dot{q}_0(t)\) scales with the dissipated kinetic energy, decreasing linearly to zero over the braking duration, while the convective term references the ambient temperature \(T_{\mathrm{amb}}\).

Quasi-static equilibrium dictates the stress field, \(\nabla \cdot \boldsymbol{\sigma} = \mathbf{0}\), neglecting inertial terms. We close the system using an orthotropic thermoelastic constitutive law. The thermal expansion tensor \(\boldsymbol{\alpha} = \operatorname{diag}(\alpha_1, \alpha_2, \alpha_3)\) is also aligned with the material axes, yielding \[\boldsymbol{\sigma} = \mathbf{C} : \bigl( \boldsymbol{\varepsilon} - \boldsymbol{\alpha}\,(T - T_{\mathrm{ref}}) \bigr), \label{eq:thermoelastic}\tag{8}\] where \(\mathbf{C}\) is the orthotropic stiffness tensor and \(T_{\mathrm{ref}}\) is the stress-free reference temperature. This reference value differs from the initial temperature if the disc contains residual thermal stress at start-up. We focus on the scalar von Mises equivalent stress, \[\sigma_{\mathrm{vm}} = \sqrt{\tfrac{3}{2}\, \mathbf{s} : \mathbf{s}}, \qquad \mathbf{s} = \boldsymbol{\sigma} - \tfrac{1}{3}\operatorname{tr}(\boldsymbol{\sigma})\, \mathbf{I}. \label{eq:vonmises}\tag{9}\] Because the stress response lags the thermal load, the two fields reach their maxima at distinct times, as defined by the peak operators in 2 . The repository accompanying this study contains the complete set of orthotropic material properties.

3.2 Finite-element reference and solver↩︎

We solve the weak Galerkin form of equations 6 to 9 on a fixed three-dimensional mesh of \(N_{\mathrm{n}} = 60914{}\) nodes. The simulation uses CalculiX [13] to perform an implicit, fully coupled temperature-displacement analysis, resolving quasi-static mechanical equilibrium via the PaStiX sparse direct solver [14]. Element definitions, time-stepping parameters, and convergence tolerances match the deposited input decks. For any input vector \(\mathbf{x}\), the solver outputs transient history fields \(\mathcal{T}(\mathbf{x},t)\) and \(\mathcal{S}(\mathbf{x},t)\), which yield the peak fields \(T(\mathbf{x})\) and \(\sigma(\mathbf{x})\) of Section 2. We compare all emulator outputs directly against this reference solution using the identical mesh layout.

The high computational cost of individual finite-element runs prevents dense, direct parametric sweeps of the input space, justifying a data-driven surrogate approach. We quantify this computational acceleration in Section 6. The simulation cost scales with the number of transient increments, which scales with the total braking duration. Higher braking loads shorten this duration, meaning simulations at high braking load require less computational time.

3.3 Design of experiments and dataset↩︎

We sample the parametric input space using a full factorial design [15]. Table 1 outlines the operational ranges and level allocations for each input parameter. This factorial grid defines the boundaries for model training and verification. While computationally demanding to construct initially, the structured grid provides uniform coverage of the parameter space and captures multi-parameter coupling. Crossing 4 loads, 4 masses, 4 initial velocities, 3 initial temperatures, and 3 ambient temperatures yields a dataset of \(576{}\) high-fidelity simulations. Generating all cases on the identical production mesh guarantees exact node-to-node correspondence.

The transient history records starting from the first saved increment, excluding the initial state at \(t = 0\). To account for variable braking durations across cases, we map the fields to a uniform time grid and mask inactive increments past the stopping time. The peak operators in 2 evaluate only this active, unmasked window.

A fixed random partition divides the \(576{}\) cases into \(403{}\) training, \(58{}\) validation, and \(115{}\) testing samples. We establish this split prior to model development and preserve it throughout the study. The test set is opened a single time, for the final metrics of Section 6, and is never used for tuning or model selection.

The dataset carries a documented representativity limit. Because the initial temperature \(T_0\) is sampled at discrete factorial levels, we use a two-sample Kolmogorov-Smirnov comparison of its empirical distribution between the test and training splits as an approximate distributional diagnostic rather than an exact hypothesis test. The p-value of \(0.0273{}\) signals a mild distributional shift, with a mean discrepancy of \(21{}\) degrees over the \(180{}\) degree range (from \(20{}\) to \(200{}\) degrees). Both splits share identical boundary limits. We treat this statistical shift as a known validation limit and analyze its impact in Section 6.

3.4 Limitations of the high-fidelity model and surrogate-to-code scope↩︎

The reference model approximates the qualified engineering setup rather than replicating the physical brake disc. Frictional heating enters the system as a prescribed boundary flux with a predetermined temporal profile, bypassing a fully coupled thermomechanical contact formulation that would dynamically link local pressure to heat generation [16]. Additionally, we assume temperature-independent material properties and neglect radiative heat transfer relative to convection. Consequently, the emulator cannot predict behavior under altered cooling environments or ventilation designs where hot-spot locations might shift [8].

These modeling assumptions restrict the physical scope of the reference. The emulator is verified directly against the finite-element solver outputs, establishing a surrogate-to-code agreement. The error budget in Section 5 measures only the distance between the surrogate and the numerical reference, leaving the physical gap to the real brake disc outside the scope of this study.

3.5 Discretization floor↩︎

The finite-element reference contains inherent numerical discretization error, representing a baseline uncertainty in the sense of verification and validation [7]. This discretization error defines a threshold below which further reduction in emulator error lacks physical or numerical meaning. We quantify this floor using three nested meshes with characteristic element sizes of \(15{}\), \(8{}\), and \(4{}\) millimeters, corresponding to \(10645{}\), \(60914{}\), and \(391480{}\) nodes. The peak von Mises stress increases from \(107.22{}\) to \(213.84{}\) and then to \(225.70{}\) megapascals, demonstrating monotonic convergence where the coarsest mesh under-resolves the stress concentrations. Because the coarse mesh lies outside the asymptotic convergence regime, we avoid Richardson extrapolation and estimate the discretization floor using only the two finer meshes. For the peak temperature, the relative difference between the production and refined meshes is \(0.002{}\), which is negligible compared to the emulator errors in Section 6. For the peak von Mises stress, the relative difference is more sensitive and case-dependent, ranging from \(0.0525{}\) to \(0.0640{}\) (equal to \(11.85{}\) megapascals at the reference case) due to localized spatial gradients. We incorporate this stress-specific discretization floor as an explicit threshold in the error budget derived in Section 5.

4 Non-intrusive POD-GPR emulator↩︎

The emulator maps parameters directly to the peak fields, bypassing solver equations and operating strictly as a data-driven model. This non-intrusive approach requires only input-output snapshots rather than full access to governing differential operators, contrasting with classical projection-based reduced-order schemes. Figure 2 outlines the two distinct phases of the formulation. During the offline phase, we construct a reduced-order basis for each physical field and train a regression model to map inputs to coordinates in this lower-dimensional space. During the online phase, we query the trained regressor at any new parameter vector \(\mathbf{x}\) to reconstruct both physical fields. We perform the spatial reduction via proper orthogonal decomposition (POD) and resolve the regression step using a Gaussian process in the coefficient space to obtain simultaneous predictions and associated variances. Although both reductions run independently, they share the same spatial discretization over the fixed finite-element mesh.

4.1 Proper orthogonal decomposition↩︎

We reduce each physical field through proper orthogonal decomposition. This process begins by centering the \(403{}\) training snapshots about their mean field and organizing them into a snapshot matrix. Singular value decomposition of this matrix yields spatial modes and associated singular values, where the squared singular values represent the modal energy [17], [18]. These modes define the directions of highest spatial variation across the training data [19]. We define the inner product using the Euclidean product of the nodal vectors, which fits our near-uniform mesh layout. Strongly graded meshes would instead require a mass-matrix-weighted inner product to account for variable element volumes. Utilizing the method of snapshots allows us to solve the underlying eigenproblem at a computational cost scaled to the training set size instead of the mesh size [17]. The resulting projection onto the first \(r\) modes yields the low-dimensional field approximation: \[u \approx \bar{u} + \Phi_r\, a, \qquad a = \Phi_r^{\top}(u - \bar{u}), \label{eq:pod-projection}\tag{10}\] where \(\bar{u}\) is the mean field, \(\Phi_r\) contains the first \(r\) spatial modes, and \(a \in \mathbb{R}^{r}\) represents the vector of reduced coordinates [20].

4.2 Rank selection by peak accuracy↩︎

Truncating the proper orthogonal decomposition basis minimizes the spatial energy norm of the reconstruction error rather than the localized peak error. These two error metrics behave differently in this setting. A low global energy error does not guarantee accurate peak predictions. The peak is a localized, low-energy feature that leading modes often fail to resolve. For instance, the first mode accounts for 0.9589 of the total energy for the von Mises field and 0.9004 for the thermal field. While the cumulative energy variance saturates using only 3 modes for the von Mises field, the localized peak reconstruction error continues to decline at higher ranks. Selecting the truncation limit based strictly on cumulative energy variance can lead to an under-resolved peak. We address this by determining the truncation rank directly from the validation set peak reconstruction error. We select the minimum rank where this peak error falls safely below our target threshold, yielding a working rank of \(6{}\) for both fields. This empirical selection remains subject to the formal verification and exact error bounds derived in Section 5.

4.3 Gaussian process regression of the reduced coordinates↩︎

We map the parametric inputs to the low-dimensional coordinates using Gaussian process regression [21]. We train an independent Gaussian process for each coordinate, requiring \(6{}\) processes per field and \(2\,6{}\) processes for the two fields in total. Treating these coordinates as independent neglects cross-modal correlations and the physical coupling between thermal and mechanical fields. This standard decoupling simplification allows efficient training, while the resulting modeling errors are captured directly by the regression term in the error budget of Section 5.

Each process maps the input vector \(\mathbf{x} \in \mathbb{R}^{5}\) to a single standardized coordinate to yield a posterior mean and a corresponding predictive variance. The covariance structure uses a Matérn kernel [22] with automatic relevance determination, assigning an independent length scale to each input parameter. We optimize these hyper-parameters by maximizing the log-marginal likelihood using multiple random restarts. Adding a minor numerical jitter to the diagonal of the covariance matrix stabilizes the linear algebra, rendering the processes near-interpolating with minimal noise. The mean field \(\bar{u}\) from 10 acts as a static offset, meaning each Gaussian process operates with a constant prior mean in the standardized space. The optimized length scales automatically scale down input dimensions that have minor influence on the coordinate predictions.

During online evaluation at a new parameter point \(\mathbf{x}\), each independent process predicts its corresponding coordinate, we reconstruct the full field via 10 , and we apply the spatial maximum operator defined in 4 . This architecture bypasses assembling or solving the finite-element system, qualifying the framework as a fully non-intrusive surrogate model [23], [24]. Alternative regression techniques include feedforward neural networks [25] or manifold interpolation methods [26]. In their standard deterministic form these do not natively output the calibrated predictive variance necessary for the a posteriori peak error budget in Section 5 and the level-one uncertainty of Section 6, although they can be equipped with uncertainty estimates by other means.

5 Peak error budget and level-one uncertainty↩︎

We evaluate the non-intrusive emulator against the high-fidelity finite-element reference at the peak state. This assessment uses an a posteriori error budget bounded by the numerical discretization floor established in Section 3.5. We also compute a level-one predictive uncertainty for the certified peak. This evaluation operates on a surrogate-to-code basis, demonstrating that the actual emulator peak error falls within the discretization uncertainty of the reference solver.

5.1 A two-term peak error budget↩︎

The certified quantity of interest is the spatial peak of each field, \(Q(u) = \max_{1 \le i \le N_{\mathrm{n}}} u_i\), defined in Section 2. Because the peak operator is nonlinear, the spatial node where this maximum occurs can vary among the high-fidelity, projected, and emulated fields. We decompose the total prediction error on this scalar quantity into two distinct sources: proper orthogonal decomposition truncation and reduced coordinate regression. Let \(Q_{\mathrm{ref}} = Q(u_{\mathrm{ref}})\) denote the peak value of the reference solver field. We define the projected peak as \(Q_{\mathrm{proj}} = Q\bigl(\bar{u} + \Phi_r \Phi_r^{\top}(u_{\mathrm{ref}} - \bar{u})\bigr)\), which represents the spatial maximum obtained by reconstructing the field from its exact projected coordinates. Let \(Q_{\mathrm{emu}} = Q(u_{\mathrm{emu}})\) be the peak value predicted by the emulator. The absolute realized error, the basis truncation error, and the coordinate regression error are formulated as: \[e_{\mathrm{real}} = \bigl\lvert Q_{\mathrm{emu}} - Q_{\mathrm{ref}} \bigr\rvert, \quad e_{\mathrm{trunc}} = \bigl\lvert Q_{\mathrm{proj}} - Q_{\mathrm{ref}} \bigr\rvert, \quad e_{\mathrm{reg}} = \bigl\lvert Q_{\mathrm{emu}} - Q_{\mathrm{proj}} \bigr\rvert. \label{eq:budget-terms}\tag{11}\] The truncation term isolates the geometric error of the reduced basis without incorporating regression inaccuracies. The regression term captures the error introduced solely by the Gaussian process models on this fixed basis. We compute these metrics using the independent test set containing \(115{}\) cases.

5.2 Error attribution↩︎

The realized prediction error on the scalar peak represents an exact algebraic sum of the signed truncation and regression contributions: \[Q_{\mathrm{emu}} - Q_{\mathrm{ref}} = \bigl(Q_{\mathrm{proj}} - Q_{\mathrm{ref}}\bigr) + \bigl(Q_{\mathrm{emu}} - Q_{\mathrm{proj}}\bigr). \label{eq:signed-decomposition}\tag{12}\] Applying the triangle inequality yields a rigorous upper bound on the absolute realized peak error: \(e_{\mathrm{real}} \le e_{\mathrm{trunc}} + e_{\mathrm{reg}}\). Case-by-case evaluation of these absolute terms across the test set confirms that this bound holds without exception. This formulation attributes the prediction error directly to its underlying mathematical sources. The basis truncation error is the primary driver of the total error, dominating the regression term in 0.86 of the thermal test cases and 1.0 of the von Mises stress cases. The linear dimensional reduction governs the prediction accuracy, while the Gaussian process approximation contributes minimally. Section 6 reports the relative peak errors.

5.3 Bounding against the discretization floor↩︎

Evaluating the realized peak error requires comparison with the baseline discretization floor of the solver established in Section 3.5 [7]. A point-by-point comparison shows that the actual emulator peak error falls below this discretization floor in 1.0 of the test cases for both physical fields. The surrogate reproduces the high-fidelity solver within its own numerical resolution limits. This alignment defines and limits our surrogate-to-code verification because attempting to achieve errors below this mesh-dependent threshold lacks physical justification.

5.4 Level-one uncertainty of the peak↩︎

We propagate the posterior coordinate variances computed by the Gaussian process regression models directly to the peak quantities of interest [21]. Using a Monte Carlo sampling procedure, we draw coordinate realizations from the posterior distributions, reconstruct the spatial fields, and extract the resulting maxima. The empirical quantiles of this generated peak distribution define a level-one predictive interval, isolating the uncertainty of the scalar peak from the full-field spatial uncertainty. This interval measures regression uncertainty alone, excluding basis truncation errors and solver discretization limits, which are handled independently within our error budget.

This propagation method exhibits two distinct mathematical behaviors. First, our Monte Carlo framework treats the coordinate posteriors as independent. The maximum off-diagonal correlation of the validation residuals reaches only 0.786 for the von Mises coordinates, and incorporating a full-covariance structure produces negligible changes in the peak coverage, verifying the validity of the diagonal assumption. Second, the maximum of a finite sample represents a downward-biased estimator of the expected peak. The predictive interval remains slightly conservative, acting as a practical engineering bound rather than a strictly unbiased statistical estimator.

We evaluate the calibration of this predictive interval on the validation split at a nominal level of 0.95. The actual spatial peak falls within the predicted interval in 1.0 of the thermal cases and 0.897 of the von Mises stress cases. The thermal predictions display slight over-coverage, while the mechanical stress predictions exhibit a minor coverage deficit. This deficit remains statistically insignificant across the validation sample because the corresponding Clopper-Pearson confidence interval, [0.788,  0.961], encompasses the nominal target. The computed level-one uncertainty remains consistent with the nominal specifications, avoiding claims of exact calibration.

6 Results on the held-out test set↩︎

6.1 Realized peak accuracy at the quantity of interest↩︎

We evaluate the emulator’s predictive performance using the spatial extremes of the temperature and von Mises stress fields over the production mesh consisting of \(N_{\mathrm{n}} = 60914{}\) nodes. To ensure a rigorous assessment, we report all subsequent metrics on an independent, held-out test set of \(115{}\) cases. This test partition remained strictly blinded until we finalized the rank-\(6{}\) emulator design, preventing any post-hoc tuning or optimization. We define a compliance threshold using a relative peak error limit of \(0.01{}\), which corresponds to the structural design tolerances of the physical brake assembly.

For the thermal peak, the emulator achieves a median relative error of \(0.000055{}\) and a maximum relative error of \(0.000438{}\) (corresponding to a maximum absolute deviation of \(0.3844{}\)), yielding a thermal compliance rate of \(1.0{}\). For the peak von Mises stress, the median relative error is \(0.00025{}\) and the maximum is \(0.017417{}\) (representing an absolute deviation of \(0.5172{}\)). This results in a mechanical compliance rate of \(0.965{}\), where only \(4{}\) out of \(115{}\) test cases exceed the target limit. We report these unadjusted compliance rates in strict accordance with our blind-testing protocol. Because the full stress field is not fully converged at a truncation rank of \(6{}\), peak identification on the truncated basis carries an epistemic risk under parametric variation. Figure 4 displays the parity plots of predicted versus reference peak values, and Figure 5 shows the cumulative distribution of the relative errors.

These performance metrics depend on two physical factors. First, the relative tolerance limit of \(0.01{}\) is highly stringent, yet the few stress cases exceeding this threshold remain significantly below the numerical solver’s discretization error described in Section 6.3. Second, because we selected the truncation rank to minimize peak error rather than global reconstruction energy, peak accuracy converges at a lower modal dimension than the full-field solution. Consequently, while the global von Mises field reconstruction does not reach a complete plateau at rank \(6{}\) (Figure 3), the emulator maintains high fidelity at the certified spatial peaks. We address the spatial localization of these extreme values in Section 7.

6.2 The two-term peak error budget↩︎

The total realized error on the scalar peak for any given test configuration equals the algebraic sum of the signed truncation and regression errors. The triangle inequality guarantees that the sum of these absolute components bounds the total absolute error, a condition verified across all \(115{}\) test cases. The finite-element discretization error is not a term of this budget. It is reported separately as a solver reference floor in Section 6.3.

In physical units, the thermal truncation error exhibits a median of \(0.0765{}\) and a \(95\textsuperscript{th}\) percentile of \(0.2726{}\), compared to a regression error median of \(0.0098{}\) and a \(95\textsuperscript{th}\) percentile of \(0.033{}\). The total realized thermal error has a median of \(0.0752{}\) and a \(95\textsuperscript{th}\) percentile of \(0.2629{}\). For the peak von Mises stress, the truncation error has a median of \(0.0197{}\) and a \(95\textsuperscript{th}\) percentile of \(0.0864{}\), while the regression error has a median of \(0.0005{}\) and a \(95\textsuperscript{th}\) percentile of \(0.0038{}\), producing a total realized stress error median of \(0.0196{}\) and a \(95\textsuperscript{th}\) percentile of \(0.0858{}\). Truncation errors exceed regression errors in \(0.86{}\) of the thermal simulations and \(1.0{}\) of the mechanical stress simulations. The spatial reduction step governs the overall emulator accuracy, suggesting that performance improvements depend on expanding the linear basis rather than optimizing the regression hyperparameters. Figure 6 displays this case-by-case error distribution.

6.3 The discretization floor and the surrogate-to-code fidelity↩︎

We evaluate the emulator accuracy against the finite-element discretization floor defined in Section 3.5 and illustrated via the mesh convergence study in Figure 7. This analysis compares a coarse mesh of \(15{}\) mm (\(10645{}\) nodes), which under-resolves the stress field with a peak von Mises value of \(107.22{}\) MPa, against our production mesh of \(8{}\) mm (\(60914{}\) nodes) and a highly refined reference mesh of \(4{}\) mm (\(391480{}\) nodes). While the thermal peaks converge fully on the production mesh, the peak von Mises stress exhibits residual mesh sensitivity. We evaluate this discretization floor, defined as the relative difference between the production and refined meshes, on a case-by-case basis. The stress discretization floor spans from \(0.0525{}\) to \(0.0640{}\). The maximum relative discretization error occurs near the median loading cases rather than the extreme loading cases because the mesh sensitivity correlates with the spatial gradient and location of the peak stress rather than its absolute magnitude. In contrast, the thermal discretization floor remains negligible, ranging from \(0.0002{}\) to \(0.0023{}\).

The total realized emulator error remains below the numerical discretization floor in \(1.0{}\) of the test cases across both physical fields. Specifically, the worst-case relative stress error of \(0.017417{}\) stays below the minimum stress discretization floor of \(0.0525{}\). When the surrogate model’s approximation error drops below the underlying discretization uncertainty of the solver, further optimization of the emulator fails to yield physically or numerically meaningful improvements. This is the basis for reading the fidelity as surrogate-to-code, under the present constant-property, radiation-free reference model and fixed mesh family. This assessment measures the emulator’s ability to replicate the numerical solver rather than the physical brake assembly, which the underlying finite-element formulation tends to overestimate.

6.4 Computational speedup↩︎

The surrogate framework replaces the high-fidelity finite-element solver (mean runtime of \(332.5{}\) seconds) with a rapid evaluation step (median runtime of \(2.238{}\) ms, excluding disk I/O for loading bases and Gaussian process hyperparameters). Since the numerical solver’s execution time scales with the physical braking duration, the computational speedup spans a wide operational range. Using a single reference machine to ensure consistent timing benchmarks, the local speedup factor ranges from \(16980{}\) to \(904399{}\), with a median acceleration of \(88921{}\). Over this benchmark suite, the total solver time of \(160271{}\) seconds is reduced by a factor of \(148579{}\). We provide the detailed distribution and range of these speedup factors in Figures 14 and 15 in 10.

6.5 Off-grid robustness and the bound on validity↩︎

We evaluate the robustness of the frozen emulator on \(15{}\) Latin-hypercube samples generated within the parameter boundaries but off the structured factorial grid. For these off-grid test points, the thermal peak prediction error exhibits a median of \(0.0210{}\) and a maximum of \(0.0744{}\), maintaining a compliance rate of \(0.2667{}\). The peak von Mises stress error shows a median of \(0.0057{}\) and a maximum of \(0.0362{}\), corresponding to a compliance rate of \(0.8000{}\). Compared to the on-grid test partition, the median thermal peak error increases from \(0.000055{}\) to \(0.0210{}\) and the compliance rate decreases from \(1.0{}\) to \(0.2667{}\). This moderate degradation reflects the challenges of off-grid extrapolation within the bounded input domain. This performance gap is more pronounced in the temperature field than in the stress field, an observation whose mechanism we do not establish here and take up in Section 8.

These metrics apply strictly within the boundaries of the sampled parameter space. This off-grid analysis serves as a secondary check to test model robustness and verify the operational limits defined in Section 3. The size of the off-grid dataset (\(15{}\) points) remains constrained by the computational cost of the high-fidelity solver rather than the surrogate, making these results indicative of model behavior under off-grid conditions.

6.6 Summary of the test-set assessment↩︎

The independent test-set evaluation yields three key conclusions. First, the emulator meets the target structural compliance criteria at the spatial peaks for the thermal field and the vast majority of the mechanical stress configurations. Second, basis truncation represents the primary contributor to the total peak error, whereas regression errors remain secondary. Third, the total prediction error remains below the solver’s inherent discretization floor across all test cases, justifying our surrogate-to-code verification framework.

7 Ablation, scope of the linear reduction, and interpretability↩︎

We analyze the primary sources of emulator approximation error, evaluate the physical limits of linear dimensional reduction, and interpret the physical dependencies captured by the surrogate model.

7.1 Ablation study and error sources↩︎

The two-term error budget in Section 5 attributes the peak prediction error to the proper orthogonal decomposition basis rather than the Gaussian process regression. Our ablation analysis confirms this attribution, particularly in the worst-performing test cases where truncation error dominates. For the minority of test cases where the reference peak von Mises stress falls outside the level-one predictive interval from Section 5.4, the regression error remains negligible compared to the truncation error. These specific cases correspond to high braking loads coupled with low aircraft masses, a regime where the localized spatial gradients are under-resolved at the selected basis rank \(r\). Consequently, the linear dimensional reduction step acts as the primary accuracy bottleneck rather than the Gaussian process regression.

7.2 Scope of the linear reduction↩︎

We assess whether adopting a nonlinear dimensional reduction scheme would improve peak prediction accuracy. Within the analyzed parameter bounds, a nonlinear reduction does not yield measurable improvements because the POD truncation error already falls below the numerical discretization floor established in Section 3.5. Any further reduction in reconstruction error would be obscured by the underlying finite-element mesh resolution, rendering more complex reduction techniques unobservable in this verification framework. The theory of the Kolmogorov \(n\)-width offers a suggestive but not conclusive backdrop [27], since it applies rigorously to linear coercive problems, whereas the present problem couples thermal diffusion with a thermoelastic response. The linear \(n\)-width of the purely diffusive thermal field decays rapidly, though we do not assume the coupled thermoelastic equations satisfy strict coercivity. Nonlinear reduction methods (including kernel principal component analysis [28], autoencoders [29], [30], and their kernel-mapped or Gaussian-process extensions [31], [32]) are designed for transport-dominated physical processes where the linear \(n\)-width decays slowly. Because our transient thermomechanical problem is not transport-dominated, we state the superfluity of a nonlinear reduction as a scoped corollary, valid for this problem class and parameter domain, and not as a general claim.

7.3 Zonal localization of the extremum↩︎

The secondary contribution predicts the extremum as a zone, the super-level set \(Z_\alpha\) of Section 2 at an operating threshold \(\alpha\), rather than as a node. We state its scope with care, since a trivial zone would carry no information. For the thermal field, this super-level set is spatially trivial at the selected threshold, spanning a complete annulus over the friction band with negligible variation across parameter cases. For the von Mises stress field, the spatial centroid and radial bounds of the zone vary less than the characteristic zone thickness, meaning the physical location of the peak zone remains too insensitive to serve as a reliable predictive output. In contrast, the total surface area of the high-stress zone represents a stable, physically informative metric. The mean area of this stress zone is \(0.3121{}\) of the total friction band surface, exhibiting a case-to-case coefficient of variation of \(0.4176{}\). This surface area displays strong monotonic correlation with aircraft mass, braking load, and initial speed across the training and validation subsets. The emulator successfully captures this physical dependency, achieving a correlation coefficient of \(0.9995{}\) on the validation set. The initial and ambient temperatures exert negligible influence on the spatial zone size. The weak correlation between zone area and initial temperature observed in the validation subset represents a finite-sample artifact present in the underlying finite-element solver data. That the emulator reproduces this local trend reflects surrogate-to-code tracking of the observed data rather than a physical causal relationship, and it remains bounded by the discretization floor, so it does not degrade the surrogate-to-code fidelity. This spatial area provides a stable scalar indicator of stress concentration that avoids the numerical instabilities associated with tracking a single maximum node.

7.4 Interpretability through the length scales↩︎

The automatic relevance determination (ARD) framework of the Gaussian process regression fits an independent correlation length scale \(\ell_d\) to each input dimension. Under the chosen Matérn kernel, these length scales provide a qualitative measure of parameter influence. Larger length scales indicate that the predicted coordinates vary slowly along the corresponding parameter dimension, signifying a lower impact on the output. Based on this analysis, the initial and ambient temperatures represent the least influential parameters, matching their minimal impact on the high-stress zone area. Conversely, the braking load, aircraft mass, and initial velocity act as the primary drivers of the physical response, which is consistent with the energy dissipation and mechanical loading mechanisms detailed in Section 3. Figure 8 displays these fitted length scales. These ARD parameters are model-dependent and differ from variance-based Sobol sensitivity indices, though both metrics yield comparable qualitative rankings [33], [34]. While Gaussian process models can be post-processed to calculate formal Sobol indices [35], and polynomial chaos expansions natively support global sensitivity analysis [36], [37], we present the ARD length scales strictly as a qualitative interpretation of parameter importance.

7.5 Section summary↩︎

Our ablation analysis establishes that proper orthogonal decomposition truncation is the dominant source of approximation error in the emulator pipeline. We verify that a linear reduced-order basis remains sufficient within the constraints of the solver’s discretization floor, and we demonstrate that the non-intrusive surrogate accurately reproduces the physical sensitivities governing the spatial extent of the stress concentrations.

8 Discussion and limitations↩︎

This section contextualizes the verification results, delineates the physical and numerical boundaries of the emulator, and positions the framework relative to existing reduced-order modeling methodologies.

8.1 What the assessment establishes↩︎

Our evaluation establishes a verified error budget anchored directly to the solver discretization floor within a formal verification and validation framework [7]. Across the independent test set, the emulator replicates the high-fidelity solver at the certified spatial peaks within the numerical discretization uncertainty of the reference model. This surrogate-to-code agreement verifies the emulator against the underlying numerical solver rather than validating it against physical experimental measurements.

Anchoring the error budget to this numerical floor establishes a rigorous stopping criterion, as optimizing the emulator beyond the spatial resolution of the reference solver yields no physical or mathematical utility. Once the emulator’s approximation error falls below this discretization uncertainty, further parameter tuning offers no actionable improvement in the fidelity of that numerical model. The two-term budget attributes this residual error. By evaluating the regression error on the fixed, truncated basis, we demonstrate that spatial truncation, rather than Gaussian process regression, governs the baseline surrogate error.

8.2 Limitations↩︎

We identify several key limitations that define the scope of these findings.

First, the peak von Mises stress reference is not fully converged on the production mesh. The stress discretization floor consequently acts as a case-dependent range rather than a static scalar, and we assert emulator compliance relative to this local interval. Although the realized peak errors on the test set consistently fall below the conservative lower bound of this range, establishing a fully mesh-independent stress reference would require resolving all training snapshots on the refined mesh, a step omitted due to computational budget constraints.

Second, the level-one uncertainty derived in Section 5.4 accounts strictly for epistemic regression uncertainty. It excludes the independent truncation and discretization errors, meaning the resulting predictive interval does not represent the full statistical uncertainty of the physical peak. The empirical stress coverage of \(0.897{}\) slightly underperforms the nominal \(0.95{}\) target, a discrepancy that remains statistically insignificant on the validation sample but could be characterized more precisely with a larger dataset. While combining these regression, truncation, and discretization uncertainties under an independence assumption would yield a single conservative interval, we do not compute such a joint metric here. This scalar budget addresses only the peak magnitude, meaning spatial localization errors represent a geometric tolerance on the stress zone area that does not affect the scalar comparisons.

Third, our assessment remains restricted to the bounded parameter space. Off-grid evaluation reveals moderate performance degradation, which is more pronounced in the thermal field than the stress field. The physical or numerical mechanism driving this asymmetric sensitivity remains unquantified, requiring a systematic sensitivity analysis of the reduced coordinates across the parameter space in future investigations. Furthermore, the surrogate’s accuracy depends on the representativeness of the offline snapshot database, and alternative operating regimes would necessitate constructing an enriched reduced basis.

Finally, the underlying finite-element formulation simplifies the physical brake physics. Neglecting radiative heat transfer, assuming temperature-independent material properties, and prescribing the frictional heat flux rather than solving a fully coupled mechanical contact problem represent major modeling simplifications. The prescribed heat flux assumption likely dominates the overestimation of the peak physical values. This work establishes code verification rather than physical validation against experimental data.

8.3 Positioning↩︎

The non-intrusive POD-GPR architecture fits within established families of projection-based and data-driven reduced-order models [5]. This framework shares structural features with non-intrusive operator inference and manifold interpolation techniques [3], [4], real-time structural surrogate modeling [11], and non-intrusive reduced-order models incorporating uncertainty quantification [38]. The primary contribution of this study lies not in the underlying surrogate architecture, but in the verification methodology that measures emulator errors directly against solver discretization uncertainties.

This methodology distinguishes our framework from adjacent reduced-order modeling paradigms. Certified reduced-basis methods provide mathematically rigorous a posteriori error bounds for field solutions or outputs [1], [2], but they require intrusive access to the solver operators, whereas our non-intrusive budget targets the scalar peak directly against a black-box solver. Bayesian calibration frameworks quantify the discrepancy between numerical predictions and physical experiments [6], addressing a surrogate-to-physics gap that remains outside the scope of this work.

Although commercial surrogate modeling toolkits offer comparable spatial field emulation capabilities, our framework explicitly accounts for the approximation error against the solver’s own discretization uncertainty. While the computational cost of building the snapshot database scales with the refinement of the high-fidelity reference model, this discretization-anchored verification protocol remains independent of mesh size. The emulator serves as an efficient computational tool, while the verified, floor-anchored error budget provides the mathematical foundation for its engineering deployment.

9 Conclusion and perspectives↩︎

The primary contribution of this work is the verification methodology itself, rather than the specific emulator architecture. It measures non-intrusive reduced-order emulator accuracy against the numerical discretization uncertainty of the reference solver. By explicitly benchmarking surrogate approximation errors against baseline mesh discretization limits, this approach extends conventional reduced-order model validation procedures. The resulting framework provides a solver-aware fidelity criterion, preventing unnecessary optimization of surrogate models below the numerical resolution already limiting the reference solution. This verification protocol remains applicable to any non-intrusive reduced-order model of parametrized finite-element simulations where a formal discretization error estimate can be established.

We demonstrated this methodology on a carbon-carbon aircraft brake disc, emulating the transient peak temperature and peak von Mises stress fields using proper orthogonal decomposition coupled with Gaussian process regression. Our a posteriori error budget successfully isolates the approximation errors associated with basis truncation from those introduced by coordinate regression. On an independent, blinded test set, the emulator replicates the peak reference solver values within the numerical discretization uncertainty of the mesh for both physical fields, confirming the surrogate-to-code fidelity of the framework. The underlying error budget reveals that linear POD spatial truncation, rather than Gaussian process regression, represents the primary bottleneck to further surrogate refinement.

Two results scope the dimensional reduction for this application. Because the linear POD truncation error falls below the verified discretization floor for the current parameter domain, adopting complex nonlinear reduction schemes, such as deep autoencoders, would provide no observable performance benefits under our verification metric [27]. This corollary is specific to the present brake family and parameter domain and does not generalize to all thermomechanical systems. As secondary contributions, the unified non-intrusive workflow predicts both physical fields simultaneously and localizes the peak von Mises stress within a probabilistic spatial zone. Although the exact coordinates of the stress maximum are highly sensitive, the surface area of this high-stress zone varies consistently with the operating parameters, serving as a stable scalar descriptor of stress concentration that the emulator accurately reproduces.

Several numerical and physical boundaries outline the scope of this study. The residual mesh dependence of the peak von Mises stress at the production resolution requires treating the stress discretization floor as a case-dependent range rather than a uniform scalar threshold. Furthermore, the level-one predictive intervals quantify only the epistemic regression uncertainty of the Gaussian processes, omitting spatial truncation and solver discretization errors which our framework handles independently. Consequently, establishing a total predictive variance would require a unified mathematical formulation combining all three error sources. The validated surrogate accuracy remains restricted to the bounded parameter space and exhibits moderate degradation under off-grid extrapolation. Finally, because the high-fidelity finite-element reference assumes temperature-independent material properties and omits radiative heat transfer, this study establishes a code-verification framework rather than a physical validation against experimental measurements.

Future research directions follow directly from these constraints. Resolving the entire offline database on the refined mesh would eliminate residual stress mesh dependencies, converting the case-by-case discretization ranges into a sharp, uniform numerical floor. Integrating the Gaussian process regression covariance, the linear POD truncation bounds, and the discretization error bands into a single, unified uncertainty framework represents an important mathematical challenge that would yield robust confidence intervals for engineering design decisions. Finally, extending this non-intrusive verification methodology from static peak predictions to full spatio-temporal histories requires addressing transient error propagation and time-dependent model reduction, determining whether our surrogate-to-code fidelity criterion remains robust in the presence of temporal error accumulation.

Data availability↩︎

Upon publication, we will make all supporting data openly available on Zenodo, including the aggregated metrics (in numbers.tex), the five-parameter design of experiments, the exact 403/58/115 data splits, and the reproducible software environment. To safeguard proprietary geometric designs, we explicitly exclude the raw four-dimensional thermal and stress tensors, the POD spatial modes, the nodal coordinates, and the underlying computational mesh. This study relies on no proprietary or client-owned datasets. Readers may contact the author to request additional details.

Declaration of competing interest↩︎

The author works within the aerospace industry. However, this study was pursued strictly as an independent, personal research project. It received no funding, resources, or involvement from the employer, who played no role in designing, executing, analyzing, or reporting this work. The author declares no other competing interests.

Declaration of generative AI and AI-assisted technologies in the writing process↩︎

The author declares that no scientific content, data analysis, mathematical formulation, or physical interpretation in this manuscript was generated by artificial intelligence. A large language model was used exclusively in a minor capacity as a linguistic aid to polish the phrasing and style of the author’s original draft. The author thoroughly reviewed all suggested refinements and assumes sole responsibility for the final publication.

10 Supplementary results↩︎

This appendix collects the figures that support the results of the main text without being essential to the argument. They are organized by theme.

10.1 Full-field reconstruction↩︎

10.2 Predictive-interval calibration↩︎

10.3 Extremum error↩︎

10.4 Computational speedup↩︎

10.5 Residual correlation and the independence assumption↩︎

10.6 Partition integrity↩︎

10.7 Zonal localization↩︎

References↩︎