Input reduction for nonlinear thermal surface loads

A multiplicity of simulations is required to optimize systems with thermal transient processes in the presence of uncertain parameters. That is why model order reduction is applied to minimize the numerical effort. The consideration of heat radiation and convection with parameter-dependent heat transfer coefficients results in a nonlinear system with many inputs as these loads are distributed over the whole surface limiting the attainable reduced dimension. Therefore, a new input reduction method is presented approximating the input matrix based on load vector snapshots using singular value decomposition. Afterward, standard reduction methods like the Krylov subspace method or balanced truncation can be applied. Compared to proper orthogonal decomposition, the number of training simulations decreases significantly and the reduced-order model provides a high accuracy within a broad parameter range. In a second step, the discrete empirical interpolation method is used to limit the evaluation of the nonlinearity to a few degrees of freedom and proper orthogonal decomposition allows the fast adaptation of the emissivity. As a result, the reduced system becomes independent of the original dimensions and the computation time is reduced drastically. This approach enables an optimal method combination depending on the number of simulations performed with the reduced model.


Introduction
The optimization of systems with thermal transient processes to reduce the energy intake or to achieve a smooth temperature distribution in order to minimize thermal stresses requires a large number of simulations. Furthermore, some system parameters like heat transfer coefficients or emissivity are uncertain. Evaluating their influence further increases the numerical effort because a wide range of parameter values has to be investigated. Against this background, a time-efficient computation of the transient temperature distribution is essential. However, in many applications the finite element method is used to discretize the thermal field resulting in large-scale models with up to several millions of degrees of freedom [1][2][3]. This leads to extensive computation times, especially if nonlinear loads like heat radiation or parameter-dependent heat transfer coefficients are considered.
For linear systems, model order reduction (MOR) is a common method to minimize the numerical costs by creating a low-dimensional model, which approximates the dynamical system properties optimally [4,5]. However, due to the nonlinearity, standard reduction approaches like the Krylov subspace method (KSM) or balanced truncation (BT), which are characterized by a high accuracy, cannot be applied directly [5]. Nonlinear reduction is the subject of current research, and the most widely used method is proper orthogonal decomposition (POD), which approximates the nonlinearity based on snapshots [3,6]. Therefore, the quality of the reduced-order model (ROM) is determined by the snapshot selection and the snapshot computation using the large-scale original model is very time-consuming [7,8]. Hence, a new method is proposed combining the advantages of KSM and POD.
The differential equation describing the discretized thermal field has the form Dθ + Cθ q(t, θ , p) (1) with the heat capacity matrix D ∈ R N ×N , the heat conductivity matrix C ∈ R N ×N , the temperature vector θ ∈ R N ×1 and the nonlinear load vectorq ∈ R N ×1 containing heat fluxes as a function of time t, temperature and an external parameter vector p [9]. The latter could, for example, express the heat transfer coefficients dependency on the velocity of the surrounding medium. N describes the number of degrees of freedom of the finite element model. The system matrices D and C are assumed to be constant.
Here, two different types of nonlinear heat fluxes are investigated: heat radiation and parameter-dependent convection. In contrast to the temperature vector, the radiative heat flux is evaluated between surface elements and not on a nodal basis. Therefore, the transformation matrices T 1 ∈ R M×N and T 2 ∈ R N ×M where M corresponds to the number of surface elements are introduced. T 1 is used to compute the vector of element temperatures θ e ∈ R M×1 θ e T 1 θ (2) and the nodal radiative heat fluxesq rad ∈ R N ×1 based on the element heat fluxesq e,rad ∈ R M×1 are calculated withq rad T 2qe,rad .
As a result of reflection, the heat flux in one enclosure depends on all surface temperatures and a system of equationsq e,rad must be solved, where σ represents the Stefan-Boltzmann constant [10]. The expression θ 4 indicates that the 4th power of the temperature is computed elementwise θ 4 i θ 4 i . The radiation matrices R 1 , R 2 ∈ R M×M with R 1 A(I − F) R 2 (I − (I − ε))F (5) depend on the view factor matrix F ∈ R M×M [10]. I ∈ R M×M is the identity matrix, while the diagonal matrices A ∈ R M×M and ε ∈ R M×M contain the element area and the element emissivity, respectively. Summarizing (2), (3) and (4), the radiative heat flux on the right side of (1) equalṡ The convective heat fluxq conv ∈ R N ×1 depends linearly on the difference between the nodal temperatures θ and the ambient temperatures θ amb that can be defined for each node individuallẏ where K ∈ R N ×N is the convection matrix [9]. As the heat transfer coefficient is adjusted for each element surface separately, the global convection matrix must be reassembled after each parameter adaptation using the modified element convection matrices. This process is much faster than the evaluation of (6) because it only includes matrix multiplications.
However, the most time-consuming step is the time-step integration of (1) and not the computation of the nonlinear loads. Therefore, the application of model order reduction to (1) is most promising in order to minimize the overall computation time. In contrast to mechanical systems, especially those used in elastic multibody simulation where model order reduction is a common method, thermal loads like convection and heat radiation are distributed over the whole surface and not restricted to a low number of interface nodes positioned at the joints [11]. As a result, thermal loads are applied to a high number of nodes and each load which can be altered independently has to be considered as a separate system input. This results in a high number of inputs limiting the achievable speedup of the computation process. Hence, the selection of an appropriate reduction method highly depends on the input number.
One of the most popular approaches for treating linear systems is KSM which performs a series expansion of the transfer function matrix and is characterized by a high accuracy and a low numerical effort as only systems of linear equations have to be solved [2,3,12,13]. Unfortunately, the reduced dimension is a multiple of the input number and consequently a high number of inputs correlates with a huge reduced dimension. Furthermore, compared to the original sparse system matrices the reduced matrices are dense and therefore the numerical effort can even rise if the number of inputs is high [14,15].
In contrast, for BT the reduced dimension can be chosen independently from the number of inputs [5,16]. This method is based on the observability and controllability matrix and those system states which are easy to observe and easy to control are retained in the reduced model [5,17]. However, for large-scale systems like those considered here the reduction process is based on the solution of the Lyapunov equations and the numerical effort herein depends on the number of inputs [4,[17][18][19].
For second-order systems like in mechanics or electrodynamics modal reduction methods are commonly used as the number of relevant modes is strictly limited due to damping that is always present in such systems [5,17]. For thermal first-order systems, the dynamic properties are not dominated by a low number of modes, especially if the local temperature distribution should be computed precisely [20,21]. Hence, a very high number of modes is necessary to achieve a good approximation and in result the modal-based computation becomes numerically ineffective.
POD, which approximates the system dynamics based on training data, is often used for the reduction of nonlinear systems [6,15,[22][23][24]. Consequently, the quality of the ROM strictly depends on the training data. A high similarity between the test simulations and the real process is required resulting in difficulties for the selection of appropriate training data especially if a large parameter range is investigated [7]. In addition, the training simulations are performed with the original model and cause an enormous numerical effort [25,26]. However, the reduced dimension does not depend on the input number and the method can be applied to nonlinear systems. That is why all publications dealing with model order reduction for thermal systems including heat radiation focus on POD [27,28].
In summary, KSM is a promising reduction method combining a high accuracy with a low numerical effort if the number of inputs is small. Therefore, for systems with a high input number a prior input reduction is required. The most commonly used approach is the singular value decomposition model order reduction (SVDMOR) which was originally developed for electronic systems and has also be applied to mechanical systems [14,29,30]. Here, similarities in the transfer functions between inputs positioned close to each other and outputs in a sufficiently large distance are used to form virtual inputs [30]. This grouping process is performed for arbitrary loads and cannot resolve local effects as it would be necessary to compute an overall temperature distribution, especially if thermal stresses should be evaluated.
Here, a new input reduction method as an adaptation of POD is presented to reduce the number of inputs instead of the system matrices allowing the application of KSM. The main advantage results from the fact that the linear parts are reduced with high accuracy and only the nonlinear load is approximated based on training data. This significantly decreases the number of training simulations and therefore the numerical effort compared to POD.
After the reduction of the temperature Eq. (1), the evaluation of the radiative heat flux (6) dominates the overall simulation time of the model because each solution step requires an expansion to the original dimension, the solution of (6) and the reduction of the load vector. That is why hyperreduction using the discrete empirical interpolation method (DEIM) is applied to restrict the computation of the nonlinear load to a low number of degrees of freedom [31]. As a result, the reduced system becomes independent of the original dimension resulting in a significant speedup.
If the emissivity is considered as an unsure system parameter, whose influence shall be evaluated based on a parameter study, the matrix R 2 in (6) which depends on the emissivity according to (5) has to be adapted for each simulation. Hence, a system of equations with the dimension M needs to be solved and the adaptation of the emissivity becomes the most time-consuming step in the simulation process. The application of POD to this system of equations allows a further reduction of the simulation time for systems with unsure emissivity.
Finally, the described reduction methods will be demonstrated using three examples for the preheating process of a pressure vessel for a large-scale geophysical experiment. Two variants use a medium sized model to investigate on the one hand the influence of the heating power and emissivity for a system considering heat radiation and on the other hand the impact of a temperature-dependent heat transfer coefficient for convection. The third example focuses on the variation of the emissivity for a large-scale model and the reduction of the computation times. The results are compared to POD.

Model order reduction
The aim of model order reduction is to approximate the original system by a much smaller one to significantly decrease the numerical effort with a minimal error in the results. Using linear MOR methods, the system is projected into a low-dimensional subspace using the projection matrix V ∈ R N ×n where n N is the degree of freedom of the reduced-order model (ROM). Applying the approach with the reduced temperature vector θ R ∈ R n×1 to (1) while multiplying with V T from the left what is called orthogonal projection leads to the reduced system [4] With the reduced heat capacity matrix D R V T DV ∈ R n×n , the reduced conductivity matrix C R V T C V ∈ R n×n and the reduced load vectorq R V Tq ∈ R n×1 (9) can be written as If n is chosen sufficiently small, the numerical effort for the time integration decreases drastically. As mentioned above, there are many different methods to compute the projection matrix V [1][2][3]5]. Here, for the reduction of the temperature vector KSM is used combining a high accuracy with moderate computation times [4,17].

Krylov subspace method
KSM approximates the transfer functions between system inputs and outputs. These inputs and output do not necessarily coincide with a single degree of freedom. Instead, they can be restricted to a selection of nodes or merge several nodes into one input with defined weight factors. That is why the load vector in (1) can be written as the product of the input matrix B ∈ R N ×L and the input vector z ∈ R L×1 q Bz (11) where L equals the number of inputs. If each node is considered as an independent input, B becomes an identity matrix. Accordingly, using the output matrix X ∈ R P×N with the number of outputs equal to P, the output vector y ∈ R P×1 can be written as y Xθ. (12) Here, the inputs and outputs of the system are assumed to be the same in order to apply an orthogonal projection and X B T . Hence, (1) transforms into the state-space formulation [4] Dθ + Cθ Bz In Laplace domain with the Laplace variable s, the outputs result as with the transfer function matrix H(s) ∈ R L×L [5]. Now, a Taylor series expansion of the transfer function matrix is performed to compute the so-called moments which shall be matched up to a predefined order between the original system and the reduced system [12]. For every expansion point, one degree of freedom is added in the reduced model for each input and each order of series expansion. Therefore, the reduced dimension n is the product of the number of inputs L, the order of series expansion o and the number of expansion points p Here, only one expansion point ( p 1) is chosen as it has proven beneficial to use a higher order of series expansion instead of multiple expansion points to achieve a specific reduced dimension because in contrast to mechanical systems the thermal system does not have a small number of dominant frequencies. The expansion point equals s 0 in order to approximate the quasistationary solution with maximum accuracy. The numerical computation of the projection matrix is based on the Arnoldi algorithm presented in [32] to project the system into a Krylov subspace because the direct evaluation of the moments is expensive and numerically unstable [12,32,33]. According to (15), the reduced dimension depends linearly on the input number. Therefore, only if L is sufficiently small, a significant reduction of the numerical effort for the simulation of the ROM can be expected. Furthermore, for the models investigated here with several tens of thousands of inputs methods like tangential interpolation for the application of the iterative rational Krylov algorithm (IRKA) have proven to be insufficient in terms of accuracy [3,13,32]. That is why the only possibility to benefit from the high accuracy of KSM is to reduce the input number prior to the application of the Arnoldi algorithm. The therefore developed approach will be presented in chapter 2.3.

Proper orthogonal decomposition
POD computes the projection matrix based on training simulations [34,35]. It is the state-of-the-art method for nonlinear model order reduction. The process is divided into an offline phase for the collection of the training data and the online phase during which the simulations of the ROM are performed [24]. As the original large-dimensional model is used for the training simulations, those comprise a high numerical effort. For the reduction of the temperature vector snapshots θ i θ (t i ) with i 1 . . . s where s is the total number of snapshots are taken at discrete simulation times t i during the offline phase. These snapshots are summarized in the snapshot matrix In the next step, the n largest singular values σ j with j 1 . . . n and the corresponding singular vectors are computed using singular value decomposition (SVD) [35] ∈ R n×n contains the singular values on the main diagonal, and the left and right singular vectors form the columns of U L , U R ∈ R N ×n . As the singular values in general decrease rapidly, n < s N can be chosen to generate a ROM with a small degree of freedom where the projection matrix equals [27] The quality and accuracy of the ROM highly depend on the appropriate choice of the snapshots and therefore the parameters used for the training simulations as well as the user experience [7,8,36,37]. As a result, many simulations of the full order model are required to guarantee the validity within a large range of input parameters corresponding to a very time-consuming offline phase.

Input reduction
As shown in (15) the reduced dimension which can be achieved using KSM depends on the number of inputs. That is why an input reduction is beneficial before the application of KSM. But it may also be necessary to reduce the numerical effort for BT. The approximation of the transfer function matrix and therefore the input-output behavior of the system requires the load vectors during simulation to be linear combinations of the columns of the input matrix B [2]. If this is matched exactly, no error occurs for the input reduction. Hence, the easiest way to fulfill this condition is to consider each degree of freedom of the model where a load is applied as an independent input leading to a large input number for distributed loads like surface or volume loads.
Loads changing synchronously are summarized to one input [38]. Then, the corresponding column of the input matrix b equals the normalized load vector Examples for such loads are heating elements where the heat flux is proportional to the applied electrical power or an equal surface temperature. The nodal heat flux at each finite element node involved changes proportional, and therefore, the ratio remains constant. More complicated load distribution can be partitioned, for example, into a mean value and an amplitude. If there are k different loads, the input matrix for synchronous loads B syn results as with j 1 . . . k. However, temperature-dependent boundary conditions like heat radiation or convection where the heat transfer coefficient is a function of temperature do not allow the representation with synchronous load vectors as each nodal heat flux changes independently. In contrast to SVDMOR, the input reduction method presented here uses training data to approximate the input matrix. The corresponding loss of generality is associated with a much higher accuracy at the load application points. The method is an adaptation of POD using load vector snapshotsq i q(t i ) with i 1 . . . s instead of snapshots of the state-space vector. Analog to (16), the snapshots are arranged in the snapshot matrix Afterward, a singular value decomposition is performed and only the l < s L largest singular values forming the main diagonal of ∈ R l×l are considered. The input matrix for the nonsynchronous loads B nonsyn equals the matrix of the left singular vectors U L ∈ R N ×l The reduced input matrix can now be used for the Krylov reduction based on (14) instead of B. Due to the drastically reduced number of inputs, a small dimension n is achieved and the computational effort for the transient simulation decreases significantly.
Only the nonlinear loads are approximated based on the training data, while the constant and synchronous loads are independent of them. Therefore, compared to POD much less snapshots are required, and the ROM is valid within a broad parameter range as it will be shown in chapter 3. Furthermore, the numerical effort for the offline simulations decreases.
However, after the reduction of the temperature equation in (10) the overall numerical effort is determined by the evaluation of the nonlinear load terms. While the computation of the radiative heat flux requires the solution of a linear system of equations in (6) for the convective heat flux, only multiplications of sparse matrices and matrix-vector multiplications are necessary. Therefore, the next chapter focuses on the reduction of heat radiation.

Reduction of heat radiation
Combining (6), (8) and (9), the determination of the reduced load vectoṙ includes the expansion of the reduced temperature vector to the original element temperature dimension M, the computation of the element radiative heat fluxes and the reduction of the load vectors. That is why hyperreduction is applied to (25) to limit the evaluation of the nonlinearity to a low number of surface elements allowing the onetime precomputation of the relevant matrices [39]. The nonlinear part f is approximated similar to (8) using the reduced nonlinear function The projection matrix for hyperreduction V HR ∈ R M×m is computed by POD, and therefore, snapshots are required [41]. As the system of Eq. (27) is over-determined, the evaluation is restricted to m rows with the Boolean mapping matrix P ∈ R M×m [15] Therefore, P contains only m elements different from zero indicating the selected degrees of freedom. If (28) is solved for f R (θ R ) and the result is inserted into (27), the nonlinear load can be approximated by Substituting f (θ R ) according to (26) in (25) with (29) the reduced radiative heat flux follows aṡ The inverse P T V HR −1 ∈ R m×m can be evaluated directly when m is chosen sufficiently small. As a result, V T T 2 R 1 R −1 2 σ εV HR P T V HR −1 ∈ R n×m can be precomputed before the simulation. Since the fourth power is calculated elementwise, the mapping matrix can be integrated in the brackets with P T T 1 V ∈ R m×n . In consequence, the terms no longer depend on the original dimensions N and M and the computation is accelerated significantly. Here, DEIM is used to determine the mapping matrix according to the algorithm presented in [31]. The selection of the surface elements considered in the hyperreduced system is iterative. At the beginning, the degree of freedom with the highest absolute value in the first column of V HR corresponding to the largest singular value is chosen. Then, during each iteration the degree of freedom leading to the largest error for the approximation of the next column of V HR is added [39]. If there is more than one enclosure for heat radiation, DEIM can be performed for each one separately and the reduced load vectors are finally summed up. Here, as the nonlinearity is evaluated elementwise the hyperreduction works straightforward. However, if DEIM is applied to convection with temperature-dependent heat transfer coefficients, the temperature of the surface element does not only depend on one nodal temperature but on the temperature of all nodes belonging to this element. That is why the number of degrees of freedom involved increases significantly. Therefore, modifications like UDEIM can be used evaluating the nonlinearity in the unassembled system [42]. As a result, the number of nodes involved and the computation time are reduced.
If the emissivity is constant, the total time for the simulation process becomes minimal. Nevertheless, considering the emissivity as an unsure system parameter whose influence on the results should be investigated requires the evaluation of V T T 2 R 1 R −1 2 σ εV HR P T V HR −1 before each simulation as the term itself depends on ε. In addition, according to (5) R 2 is a function of the emissivity. Hence, a new factorization of R 2 is necessary after each adaptation limiting the speedup. That is why a third reduction is applied to the system of equations with the corresponding solution vector x ∈ R M×1 . As the system is only solved once for each emissivity, classical reduction methods for linear systems like Guyan reduction or KSM do not decrease the simulation time because the computation of the projection matrix would take longer than solving the system of equations once. Therefore, POD is used with the projection matrix W ∈ R M×m with m ≈ m M and the reduced form of (32) becomes Inserting (33) in (30), the radiative heat flux finally can be computed aṡ Now, the transient simulation and the adaptation of the emissivity can be performed at minimal numerical costs.

Numerical examples
Three examples are used to demonstrate the advantages of the described input reduction method in combination with KSM. All of them focus on the optimization and parameter variation for the pressure vessel of the DRESDYN (DREsden Sodium facility for DYNamo and thermo-hydraulic studies) project [43][44][45][46]. The aim of this large-scale experiment at the Helmholtz-Zentrum Dresden-Rossendorf is to investigate the development of planetary magnetic fields due to rotation and precession. The central unit of the experimental plant, which is shown in Fig. 1, is a pressure vessel with an inner diameter of 2 m rotating around its longitudinal axis with up to 10 Hz. Furthermore, the platform on which the vessel is installed rotates around a second axis orthogonal to the ground with a maximum frequency of 1 Hz. The resulting gyroscopic effects drive the flow of the liquid sodium within the vessel mimicking the conditions in the outer core of the earth for example. As sodium has a melting temperature of 97°C, the experiments will operate in the temperature range between 110 and 170°C.
The sodium is stored and melted in external tanks. That is why the vessel must be preheated prior to the filling process to avoid solidification which could cause an imbalance because of gas locks. Therefore, all surfaces being in contact with the sodium should have a temperature of at least 110°C. The preheating is done  via heating foils, which can only be applied to surfaces not in contact with the sodium. The heating power can be regulated individually for three sections (q 1 ,q 2 andq 3 in Fig. 2). The volumes 1, 2 and 4 will be filled with sodium during operation, while volume 3 contains argon in a compressible bellows to compensate the thermal expansion of the liquid sodium. However, during the preheating process all four volumes are filled with argon to avoid chemical reactions between the later filled in sodium and the air moisture. Volume 5 is located behind the sealing and contains air. The vessel does not rotate during the preheating, and heat transport due to convective flows in the gas can be neglected. Due to the restricted positioning of the heating foils and the low thermal conductivity of the stainless steel used for the vessel's hull, the heat supply, especially for the interior structures, is dominated by radiation. Therefore, the consideration of heat radiation in the five described enclosures becomes essential. Furthermore, at all outer surfaces heat radiation and convection to the environment are considered and at the bearing, a temperature boundary condition is defined. Unless otherwise stated, the load parameters summarized in Table 1 are used. All other parameters vary for the three applications. The transient simulation is always performed for a time interval of 10,000 s. Due to symmetry, the model is limited to the sixteenth part of the pressure vessel. For all simulations, the maximum error should be less than 1 K.

Variation of heating power and emissivity
The aim of the first example is to evaluate the validity of the ROM within a broad parameter range at a minimum number of training simulations. That is why the heating powers and the emissivity are varied within the intervals described in Table 2. Therefore, the 100 parameter combinations illustrated in Fig. 3 were generated randomly. All other parameters remain constant. As a result, the convection matrix can be integrated in the conduction matrix.
The number of training simulations is increased iteratively. In the first iteration, only one transient simulation (q 1 3 W/cm 3 ,q 2 0.5 W/cm 3 ,q 3 2.5 W/cm 3 , ε 0.5) marked by the larger light green circle in Fig. 3 is used to collect the snapshots. Based on these data, model order reduction is performed and the error for each of the 100 parameter combinations is evaluated in comparison with the original model. The parameter combination resulting in the largest error is added to the training data for the next iteration.   Table 3 summarizes the dimensions of the original and the reduced-order models. The largest of the five enclosures contains 5125 surface elements radiating toward each other, while M ∞ describes the number of elements with radiation to the environment. Considering the three different heating powers, the temperature boundary condition and the convective heat flux from the environment there are k 5 synchronous loads. For the reduction of the radiative heat flux l 20, singular vectors are taken into account leading to a reduced dimension n 200 according to (15) since the order of series expansion is chosen as o 8. For the reduction of the heat radiation using DEIM and POD as well as the reduction of the temperature equation with POD, all singular values larger than the computation accuracy of the system are retained with a maximum of 200. Hence, the reduced dimensions increase with each iteration as more snapshots are considered. Therefore, Table  3 indicates the upper limits.
Combining the new input reduction method with the Krylov reduction (KSM in Fig. 4), the maximum error is already below the intended limit for the first iteration. Therefore, one training simulation is sufficient to approximate the whole parameter range. In comparison, using POD six iterations are required and especially at the beginning enormous errors occur showing that it is necessary to sample the parameter space carefully. For both methods, the computational effort for the time integration is comparable and dominated by the evaluation of the heat radiation. If in addition to the Krylov reduction the heat radiation equation is also reduced applying DEIM (KSM + DEIM) or DEIM and POD (KSM + DEIM + POD), the computation time can be minimized. However, more training simulations are required to achieve the desired accuracy, but in comparison with POD is considered as a function of temperature, while heat radiation is neglected (ε 0). One training simulation is performed (q 1,train 3 W/cm 3 ,q 2,train 0.5 W/cm 3 ,q 3,train 2.5 W/cm 3 ), and afterward, the heating powers are varied simultaneously with the load factor ξ As the evaluation of the nonlinear load requires only matrix multiplications, no hyperreduction is applied in this example. The dimensions are summarized in Table 4. Here, only the three heating powers and the temperature boundary condition can be considered as synchronous loads (k 4) and during the input reduction the singular vectors corresponding to the l 13 largest singular values remain in the input matrix. With an order of series expansion o 8 in the Krylov reduction, the reduced dimension becomes n 136 according to (15). For POD, all singular values smaller than the floating-point accuracy are truncated.
Despite a large load variation, only minimal errors occur using KSM, while with POD the error is two orders of magnitude larger (Fig. 5). Therefore, small deviations of the load lead to significant errors and more training simulations would be necessary to achieve the required accuracy.

Large-scale model
To limit the number of simulations using the original model, here only the emissivity is varied in the interval ε 0, 0.05 . . . 1 (q 1 3 W/cm 3 ,q 2 0.5 W/cm 3 ,q 3 2.5 W/cm 3 ). The model uses a finer mesh in order to compute the temperature distribution more precisely resulting in a degree of freedom which is significantly larger compared to the other examples (Tables 5 and 6). Analog to chapter 3.1, five synchronous loads are considered. In addition, l 11 singular values remain in the reduced input matrix resulting in a reduced dimension n 128 for the Krylov reduction as the order of series expansion is chosen to o 8. Two transient simulations ε 0.2 and ε 0.8 are performed to compute the snapshots.
Similar to the other examples, POD corresponds to the largest errors, especially for ε 0 (Fig. 6). In contrast, KSM produces very accurate results if the heat radiation is neglected because the linear loads are approximated without snapshots. For POD, at least three training simulations are necessary to reduce the maximum error below 1 K, while KSM with unreduced heat radiation requires only one training simulation. If in addition to the temperature equation also the heat radiation equation should be reduced using DEIM or DEIM and POD, the error increases slightly and with two training simulations the desired accuracy is achieved. Here, significantly lower deviations between the original model and the ROM can be observed when only DEIM is applied.   Table 7 compares the numerical effort for the adaptation of the radiation matrix using (32) or (33) and for the time integration. Reducing the temperature equation, the simulation can be accelerated by factor 40. DEIM allows a further decrease of the simulation time to only 0.2 s. If the emissivity remains constant and the radiation matrix must not be adapted before each simulation, the computational effort during the online phase is reduced by four orders of magnitude. Otherwise, the time required to solve (32) becomes dominant and applying POD according to (33) minimizes the time for one simulation cycle to approximately 3 s instead of 1 h. The additional computational effort during the offline phase to collect the required snapshots only proves advantageous if the simulation must provide real-time data or if the total number of simulations is sufficiently high. For the example shown here where the reduction of the heat radiation requires two training simulations instead of one, at least 42 different variants have to be investigated combining KSM and DEIM and 28 when POD is used in addition for the adaptation of the radiation matrix. POD for the temperature reduction delivers no significant advantage in terms of simulation time. However, this method can also be combined with the reduction of the heat radiation. Then, the error is dominated by the temperature equation and the accuracy is similar to the model with unreduced heat radiation.

Conclusion
The consideration of heat radiation or temperature-dependent heat transfer coefficients for convection leads to a system of nonlinear equations for the temperature field inhibiting the application of standard model order methods for linear systems. That is why the nonlinear terms are treated as an external load instead of integrating them into the conduction matrix. As a result, the system matrices become constant. However, besides the temperature equation a second system of equations must be evaluated to compute the load vector resulting from the current system state. Therefore, the problem investigated here can be interpreted as a coupled field analysis using load vector coupling. On the one hand, the temperature field determines the right-hand side for radiation equation, while, on the other hand, the radiative heat flux is a part of the thermal load vector.
Both fields interact at many nodes, in this case at all surface nodes with applied convection or heat radiation. Hence, the number of system inputs becomes large limiting the application of some model order reduction methods. As the solution of the temperature equation has proven to be most time-consuming during the transient simulation due to its high dimension, in the first step only this system was reduced. A new input reduction method was presented approximating the input matrix based on the singular value decomposition of load vector snapshots. The reduced input matrix can then be used in combination with standard reduction methods based on the state-space formulation like KSM or BT.
The advantages of the new input reduction method combined with KSM were demonstrated for three numerical examples. Compared to POD, which is the state-of-the-art method for such problems, the number of snapshots necessary to generate a valid ROM was drastically reduced. Therefore, the numerical costs during the offline phase of the simulation are minimized. Furthermore, despite large differences between the training data and the transient analysis the error between the original model and the ROM did not grow extensively. That is why the new method provides more reliable models especially if a broad parameter range is investigated.
In a second step, DEIM was used to limit the evaluation of the nonlinearity, which in case of heat radiation results from the fourth power of the temperature, to a small number of degrees of freedom. As a result, the numerical effort during the transient analysis decreased by four orders of magnitude for a large-scale model. Finally, POD was applied to the radiation equation to accelerate the adaptation of the radiation matrix at changing emissivity. The reduction of heat radiation requires more training simulations than the input reduction to achieve the desired accuracy, but compared to POD significantly less snapshots are necessary. Therefore, the three steps allow the customization of the computation process depending on the number of simulations using the ROM performing parameter variation, optimization or uncertainty quantification at minimal numerical costs.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest No potential competing interest was reported by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.