Reduced order homogenization of thermoelastic materials with strong temperature dependence and comparison to a machine-learned model

In this work, an approach for strongly temperature-dependent thermoelastic homogenization is presented. It is based on computational homogenization paired with reduced order models (ROMs) that allow for full temperature dependence of material parameters in all phases. In order to keep the model accurate and computationally efficient at the same time, we suggest the use of different ROMs at few discrete temperatures. Then, for intermediate temperatures, we derive an energy optimal basis emerging from the available ones. The resulting reduced homogenization problem can be solved in real time. Unlike classical homogenization where only the effective behavior, i.e., the effective stiffness and the effective thermal expansion, of the microscopic reference volume element are of interest, our ROM delivers also accurate full-field reconstructions of all mechanical fields within the microstructure. We show that the proposed method referred to as optimal field interpolation is computationally as efficient as simplistic linear interpolation. However, our method yields an accuracy that matches direct numerical simulation in many cases, i.e., very accurate real-time predictions are achieved. Additionally, we propose a greedy sampling procedure yielding a minimal number of direct numerical simulations as inputs (two to six discrete temperatures are used over a range of around 1000 K). Further, we pick up a black box machine-learned model as an alternative route and show its limitations in view of the limited amount of training data. Using our new method to generate an abundance of data, we demonstrate that a highly accurate tabular interpolator can be gained easily.


Introduction
The consideration of anisothermal settings is mandatory during the simulation of composite materials generated via laser melt injection [10,24]. Such composite materials may be obtained by adding fused tungsten carbide particles (FTC), for instance, into the melt pool of a metallic substrate to generate surface coatings [39]. These coatings are used to boost wear resistance, more precisely to protect metallic surfaces against abrasion and erosion, e.g., Zou et al. [39]. They can significantly extend the lifetime of engineering parts due to the outstanding material characteristics of the locally produced metal matrix composite (MMC) while retaining the ductility of the substrate.
The thermomechanical response of such microheterogeneous materials has to be investigated to understand the response on the macroscopic scale, i.e., on the scale of components or parts. Obvious influence factors on the microscopic response are, for instance, the phase volume fraction and the shape and orientation of the inclusion S. Sharba · J. Herb · F. Fritzen (B) Data Analytics in Engineering, Institute of Applied Mechanics (MIB), University of Stuttgart, Universitätsstraße 32, 70569 Stuttgart, Germany E-mail: fritzen@mib.uni-stuttgart.de phase. Various methods were developed to predict the linear thermoelastic response of microstructures, e.g., the seminal work of Eshelby [6] that investigates the special case of an ellipsoidal inclusion and its extension by Mori and Tanaka [30] to predict the response of structures with near isotropic inclusions. Other methods aiming at finding rigorous bounds of the effective properties of composites were developed by Hashin and Shtrikman [18] and generalized by Willis [38]. While the aforementioned approaches tackle elastic problems, semianalytical methods allow the estimation of the effective behavior of some nonlinear and inelastic composites [3,23].
Besides analytical and semi-analytical approaches, computational homogenization is commonly employed. It can accurately predict the thermomechanical response of heterogeneous structures, e.g., Feyel [8], Feyel and Chaboche [9] and Miehe [29]. In such schemes, the constitutive response of every material point on the macroscale is obtained by solving a microscopic boundary value problem (BVP) on a representative volume element (RVE) with boundary conditions from the corresponding macroscopic point. Direct numerical simulation (DNS) on the RVE scale, e.g., in the FE 2 context, is computationally demanding and usually prohibitive due to its long computational time and large memory demands as well as due to its massive energy footprint. On the contrary, reduced order models (ROM) and data-driven surrogate models provide an appealing and efficient alternative to DNS. The goal of ROMs is to provide an accurate and efficient computational scheme by reducing the large number of degrees of freedom on the microscopic scale systematically, i.e., in a problem-specific way. Reduced order computational homogenization approaches to predict the linear response of composites are well-established by now, e.g., Fritzen et al. [14] and Fritzen and Böhlke [11]. However, when it comes to temperature-dependent microscopic material behavior, the task gets more challenging and nontrivial. For each temperature observed in a macroscopic simulation, the underlying microscopic constituents within the composite show a different thermomechanical response.
In the tackled problem of laser melt injection, temperature changes cannot be ignored due to the very wide temperature range of the process, i.e., from room temperature (around 300 K) to melt temperature (around 1356 K) within seconds (and vice versa). This rapid variation of the process temperature strongly affects all thermomechanical material parameters.
Hence, our goal in this paper is to develop a reduced order model to tackle the challenging anisothermal setting. The developed approach is built on field solutions (displacements, strains, stresses) obtained via DNS at (few) discrete temperatures. Starting from this sequence of ROMs-each of them being exact at their respective temperature-we introduce a new but very small optimization problem at arbitrary intermediate temperatures.
The latter can be solved in real-time yielding accurate homogenized properties. Further, the proposed method allows the reconstruction of the full local displacement, strain, and stress fields in an inexpensive postprocessing step which enables statistical failure analysis and more.
Regarding DNSs, an in-house FE simulation software called Combo FANS was used. Therein, a Fourier-Accelerated Nodal Solver (FANS) [25] is used to iteratively solve the linear system of equations by successive discrete convolutions of residual forces with fundamental solutions matching discrete Green's functions. The fundamental solutions are available in closed form in the Fourier domain. Further, the convolutions can efficiently be computed in the frequency domain before a computationally sleek inverse fast Fourier transform (FFT) is used to update the unknowns. In addition and for the sake of efficiency, Combo FANS utilizes a coarse graining technique based on composite voxels [15,21,22]. Inside of coarse-grained voxels comprising a material interface, the material behavior is replaced by that of an infinite laminate with matching orientation and phase volume fraction. Recently, the usefulness of advanced 3D image processing and a generalization of composite voxels toward non-equiaxed so-called composite boxels (ComBo) have been illustrated in Keshav et al. [22]: The image processing yields more accurate orientation information, while ComBo enables axisspecific coarse graining which can lead to improved computational performance for pseudo-uniaxial materials, containing almost aligned fibers. Besides a much improved effective material response, it is noteworthy that the local stress fields benefit significantly. Stair-casing is reduced, and accurate information on interface tractions can be extracted.
In view of reduced order models, we focus on classical Galerkin projection which is established today, see, e.g., the textbooks [19,[31][32][33]. It is well-known that the linear homogenization problem can be solved by the superposition principle for arbitrarily applied strains. However, such ROM, exact at its identification temperature, has limited validity for different material data. Naive interpolation assuming a convex combination of the solutions at different temperatures would be the go-to engineering approach. However, if the number of sampling temperatures is small (e.g., in order to limit the amount of data as well as to have low offline cost), such naive interpolation can lead to inaccuracies. Hence, we suggest a recombination of the available basis functions at select, i.e., corresponding to a few temperatures, to yield an energy optimal basis that is local with respect to the temperature.
The outline of the paper is as follows: In Sect. 3, the homogenization problem is introduced. The proposed reduced order homogenization methods are described in Sect. 4, and algorithmic details are summarized in Sect. 5. The numerical results for different benchmark problems are presented in Sect. 6. Section 7 provides a concise summary and concluding remarks.

Notation
An index-free notation is used for different tensor types. Scalars are represented by lowercase letters (a, ϕ). First-order tensors are denoted by lowercase boldface letters, e.g., (u, x) while uppercase or Greek boldface letters (I, N, σ , ε) refer to second-order tensors. Blackboard bold letters denote fourth-order tensors (C, H). Discrete values are stored in column vectors (•) or matrices (•). Symmetric second-order tensors such as stress σ or strain ε are identified by their corresponding six-dimensional column vector, in Mandel notation which is detailed in section A.. The volume average of a general local field on a representative volume element (RVE) is defined by where | | denotes the volume of the domain . If a function f is dependent on a strain field ε and it is evaluated at a given temperature θ , then it is expressed via f (ε θ), i.e., θ is considered an external parameter.

Homogenization problem
The proposed approach is introduced within a first-order homogenization scheme of microheterogeneous materials with anisothermal constitutive relations in an infinitesimal strain framework. In view of mechanical homogenization, separation of length scales is assumed, i.e., gradients of macroscopic fluctuations are assumed to be small in comparison with the characteristic length of the representative volume element (RVE). The macroscopic domain is denoted by , and its corresponding macroscopic fields are overlined (•). Balance of linear and angular momentum must hold on either of the two scales: where σ and σ are the stress fields, b denotes a volumetric force density and BC and BC refer to boundary conditions on the macroscopic scale and on the RVE, respectively. The scale coupling relations are defined via the volume averaging operator (assuming no voids are present) as The boundary conditions of the microscopic boundary value problem in Eq. (3) are taken in accordance with the Hill-Mandel condition [20]. The Hill-Mandel condition is, for instance, satisfied by periodic displacement fluctuations u, i.e., the microscopic displacements u, and the strain field ε can be written in terms of their macroscopic counterparts as where x ∈ is the position in the microscopic domain, and the infinitesimal strain tensors ε and ε are defined via ε = u ⊗ sym ∇, ε = u ⊗ sym ∇. The superimposed rigid body translation is represented by u, and infinitesimal macroscopic rotations are eliminated, therefore, the use of ε. The periodic displacement and strain fluctuation fields are denoted by u and ε, respectively. To ensure the clarity and completeness of this work, important definitions are listed below.
Definition 2 A statically admissible or self-equilibrated stress field σ is a field that satisfies the balance of linear and angular momentum in the absence of body forces (b = 0), i.e., it is a divergence-free field as in eq. (3).
Definition 3 A compatible strain field or kinematically admissible strain field ε = u ⊗ sym ∇ is a field equal to the symmetric gradient of a displacement field u that satisfies the kinematic boundary conditions.

Definition 4
Strain fluctuation is the symmetric gradient of a displacement fluctuation field, i.e., it is a compatible field in itself. Further, the volume average of such a field vanishes and does not affect the corresponding strain average. Hence, any multiple of such a field can be added to any kinematically admissible field without introducing any incompatibilities and without changing the applied strain ε.
Definition 5 A periodic displacement-fluctuation field is a field that has identical displacement values on opposing boundaries of an RVE, and the volume average of such field is zero.

Definition 6
The decomposition of microscopic strain states that given an arbitrary strain ε, it can be additively decomposed into a homogeneous part of a virtually homogeneous material ε and a strain fluctuation part ε via ε = ε + ε.
The microscopic temperature θ is written analogously to the displacement but asymptotically, i.e., as the RVE length tends to zero l → 0, it coincides with θ . Hence, temperature-dependent material parameters in the microscopic scale are considered as functions of the macroscopic temperature θ , e.g., Chatzigeorgiou et al [4], which is in line with the assumed separation of length scales. Therefore, when seeking the mechanical response in the sequel, the temperature θ is assumed to be constant across , i.e., θ = θ which can be seen as an external parameter.
In this work, the focus is confined to anisothermal thermoelasticity, where material parameters are temperature-dependent. Therefore, on the microscopic scale, the anisothermal elastic material response is assumed to be driven by a Helmholtz free energy density function where C(θ ) is a symmetric and positive definite material stiffness tensor, ε θ (θ ) is the thermal strain at temperature θ and ψ θ (θ ) is the thermal energy. The stress σ (ε θ) and the stiffness C(θ ) are derived via the relations The thermal strain is defined as the integral of the thermal expansion coefficient α(θ ) through where θ 0 is the reference temperature at which the thermal strain vanishes, i.e., it refers to a stress-free ground state with ε θ (θ 0 ) = 0.

Solution of the thermoelastic homogenization problem
In the homogenization problem introduced earlier in this section (cf. Equation (3)), the microscopic problem has to be solved for every material point on the macroscopic scale and for every given loading condition (ε, θ) in order to obtain σ (or also C) and, possibly, locally resolved stress and strain fields which could be used, e.g., for failure initiation indication. Hence, a true many-query context is present which would rule out direct numerical simulations due to the related computational cost.
In the following, we reformulate the microscopic problem that has to be solved as where p is a generic eigenstrain (or polarization) that, in the case of thermoelasticity, comprises the macroscopic strain ε (constant throughout ) and the heterogeneous, temperature-dependent thermal strain ε θ (θ ) according to In the following, the vector and tensor notation are used interchangeably, e.g., p ↔ p. For stress and strain tensors, the orthonormal basis described in section A is assumed. Then, we can rewrite Eq. (11) through The matrix P(θ ) contains the macroscopic and thermal eigenstrains, while the vector ζ represents the activation coefficients for these modes. To this end, only the thermal eigenstrain ε θ is variable. It depends on the temperature and the material at position x.
The general formulation (10) with the split of the total strain ε into a fluctuation part ε and a polarization or eigenstrain p simplifies the following derivation. It also renders the derivation more general, e.g., for future applications to elastoplasticity where the eigenstrain in Eq. (11) will have additional contributions emerging from the plastic strain.
Effective RVE properties, which are required for the overarching macroscopic simulation, result from solutions to the microscopic problem. By virtue of the linearity of (10), these solutions can be obtained through the superposition principle: At a given temperature θ = θ , Eq. (10) is solved 7 times to obtain displacement fluctuation fields u 1 , . . . , u 6 , u θ and related strains ε 1 , . . . , ε 6 , ε θ . More precisely, Eq. (10) is solved for ζ β ← e β (β = 1, . . . , 7). This corresponds to solving (in the absence of thermal eigenstrains) the elastic homogenization problem with temperature-dependent elastic parameters and, further, the thermal eigenproblem in the absence of macroscopic deformation. The resulting thermoelastic displacement, strain, and stress are obtained from: The operators U , E and S are given by The strain and stress localization operators E, S contain compatible strains and admissible stresses (cf. Definitions 3 and 2). They can be written as Further, the volume average over the RVE has the following properties: where τ θ (θ ) is the effective stress from the deformation constraint thermal eigenproblem.

Effective material model
The material belongs to the class of generalized standard materials (GSM) as in Suquet et al. [36]and Halphen and Nguyen [17]. By using a Galerkin formulation, the weak form of Eq. (10) can be obtained from the free energy as a function of the strain fluctuation ε Minimizing this energy with respect to the strain fluctuation ε results in the stationary condition for where δ ε denotes kinematically admissible variations of ε. The stationary point of (i.e., its minimum) is the effective strain energy ψ Defining the strain operator by B ε and the discrete degrees of freedom (nodal DOF or reduced kinematic DOF) by u, the relation is gained. Condition (22) leads to the following algebraic equation: Note that in Eq. (25) and thereafter dependencies on x, ε and θ are not written explicitly to simplify the notation. In general, i.e., in a nonlinear setting, the system in Eq. (25) can be solved by a Newton-Raphson scheme with the following Jacobian: More precisely, J coincides with the FE stiffness matrix in a finite element context and with its reduced counterpart in a Galerkin reduced order model (ROM). The effective stress and tangent stiffness are deduced from Following established procedure (see, e.g., [12]), the effective stiffness is given by The symmetry of C is guaranteed by construction, and positive definiteness can be confirmed as in Fritzen and Kunc [13]. The first term in Eq. (28) is the upper Voigt/Taylor bound that overestimates the tangent stiffness, while the second term accounts for the fluctuations and leads to a relaxation of the overly stiff estimate. Notably, this part depends on the admissible strain space, and it is valid for finite element (FE) models as well as ROM, where loose inequalities relate the spectrum of C (ROM is stiffer or equal than FE), since the only difference is found in the operator B ε .
At either of the two discrete temperature indexed i ∈ {1, 2}, the procedure from Sect. 3 yields operators The superscript is also used to indicate other fields, e.g., the temperature-dependent stiffness . Balance of momentum Eq. (10) is satisfied at both temperatures for arbitrary applied loadings ζ = [ε, 1], i.e., At the intermediate temperature θ (α) , the stiffness C (α) and the thermal strain ε (α) θ are defined by the underlying constitutive models. The target of this paper is the construction of near exact, yet approximate solution operators U (α) , E (α) , S (α) that depend only on the available solutions at the boundary of the considered interval, namely θ (1) and θ (2) . Naive linear affine combinations of E (i) or S (i) are, however, prone to fail: • While affine combinations of the strain fields preserve kinematic admissibility, the temperature-induced variation of the stiffness generally induces stresses that do not satisfy the balance of linear momentum Eq. (10). • Likewise, affine combinations of stresses preserve static admissibility, but the thermal variation of C (α) would-in general-lead to incompatible strains ε incomp = C −1 σ + ε θ .

Proposition
As stated before, it is nontrivial to solve Eq. (10) in closed form due to the variability of both C (α) and ε (α) θ . While it is feasible to solve Eq. (10) at specific discrete temperatures given a certain computational investment. Finding a generic solution for arbitrary intermediate temperatures θ (α) with θ (1) ≤ θ (α) ≤ θ (2) remains challenging, even when the solutions at θ (1) and θ (2) are known. In the sequel, we propose a thermoelastic homogenization procedure that makes use of (few) solutions at discrete sampling temperatures. These solutions serve as samples/snapshots that are used to predict solutions at arbitrary intermediate temperatures in a real-time fashion. The task is thus to get optimal operators U (α) , E (α) , and the induced S (α) in an efficient way through (non-Naive) recombination of the available data, see Fig. 1.
Since U and E are linearly related and S is an inferred quantity, we focus on finding E (α) . In order to a priori guarantee kinematic compatibility (see Definition 3), we propose a linear combination of the strain Fig. 1 Visualization of the inputs and outputs of the proposed interpolation method. The effective total strain is shown as a representation of E fluctuations at either temperature which can be defined by matrices φ (1) Thus, we seek a matrix φ (α) ∈ R 14×7 that provides a temperature-dependent, optimal interpolation of available field data. The stress localization is obtained using Eq. (18) via Note that φ (α) depends on the temperature, but not on the spatial coordinate x or on the applied strain ε, while E (i) (i = 1, 2) is a discrete spatially heterogeneous operator at θ (i) . Finding φ (α) is the core of the proposed algorithm. Based on Eqs. (31) and (32), displacement, strain, and stress fields are obtained via where (C (α) , P (α) , φ (α) ) are the only temperature-dependent quantities. The missing interpolation matrix φ (α) is gained from the minimization problem (23). First, we express the potential Eq. (21) at θ = θ (α) as The stationary condition for gaining the effective free energy (23) results in the matrix-valued root-finding problem After closer inspection, this algebraic system reduces to Allowing ε to be arbitrary and accounting for the structure of ζ cf. (12), we find that in order for (38) to hold, ζ can be ruled out. Notably, this system could have a nonunique solution since K (α) ∈ Sym ≥0 (R 14 ), i.e., K (α) can be positive semi-definite. Therefore, the solution is obtained via least-squares solution or, equivalently (up to numerical issues), by using the Moore-Penrose pseudoinverse K (α) † of K (α) : The solution to the problem thus consists of the following steps: • Precomputation of solutions at (few) discrete temperatures (offline) • Assembly of K (α) ∈ R 14×14 and F (α) ∈ R 14×7 (online) • Solution of the least squares problem (39) (online) • Postprocessing of the solution (online) (e.g., effective stiffness/thermal expansion, local stress fields, …).
Here, offline refers to preprocessing calculations effected prior to the evaluation of the model and online refers to operations required during the model query.

Stiffness and thermal strain approximations
In Sect. 4.1, it is assumed that the material stiffness C(x | θ) and the thermal strain ε θ (x | θ) are exact. Hence, the assembly of the matrices K (α) and F (α) which depend on the temperature θ (α) during the online phase is the computational bottlenecks of the method. In this section, we refine the scheme by making use of an approximate affine parameter decomposition. The temperature-dependent coefficients of the problem emerge from the variation of the stiffness C where SI X (m) (β) is a short-hand notation for the secant interpolation of quantity X in phase m at secant parameter β. Hence, the 2M free parameters of the secant interpolation are the secant interpolation coefficient for the stiffness α C (m) and for the thermal strain α ε θ (m) within each phase m as a function of given α θ (see (29)). In the sequel, we investigate four different approaches which differ in how these coefficients are chosen: • Approach (O 0 ): One explicit DOF (α θ ) and no implicit DOF The coefficients match the one of the temperature, i.e., α C • Approach (O 1 ): One explicit DOF (α θ ) and one implicit DOF One new coefficient α is introduced, i.e., α C (m) = α ε θ (m) = α: • Approach (O 2 ): One explicit DOF (α θ ) and two implicit DOFs One coefficient for the stiffness and the thermal strain of all materials are introduced (α C , α ε θ ): • Approach (O 2M ): One explicit DOF (α θ ) and 2M implicit DOFs All coefficients in all phases are defined independently: The coefficients α (O 1 ), α C , α ε θ (O 2 ) and α C  Table 1 for a microstructure with two constituents, i.e., M = 2. O 2 is recovered from O 2M ↔ O 4 by dropping the material number subscripts, while O 1 is recovered by dropping all subscripts and superscripts.
In addition to using up to 2M values for secant interpolation of the material data to compute φ (α) , an alternative naive interpolation can be studied. It corresponds to an approximation of the strain localization operator via In the sequel, this is labeled (N). It corresponds to settings φ (α) = (1 − α θ )I , α θ I T , i.e., a predefined (generally not optimal) choice for the interpolation of the reduced basis of the strain is made (see also the comment at the of Sect. 4.1).

Affine operator representation and efficient implementation
The implementation is detailed for the most general case, i.e., O 2M with 2M parameters. As detailed above, the O 2M approach contains the other cases. For an intermediate temperature θ (α) ∈ [θ (1) , θ (2) ] characterized by α θ ∈ [0, 1], the eigenstrain from Eq. (11) is evaluated with the help of Eq. (41) as where • = • (2) − • (1) . Likewise, we reformulate (40) as . (48) The linear system Eq. (39) is rewritten as where all terms can be precomputed given quantities corresponding to the temperatures θ (1) and θ (2) at the boundary of the considered interval via Given these preassembled matrices, Eq. (49) provides the condensed stiffness as a linear function of the interpolation parameters, while the right-hand side F (α) is bilinear in the coefficients. The assembly and the solution of the least squares problem are feasible in real time.
Likewise, an approximate stress localization operator from Eq. (32) can efficiently be assembled via The last line emphasizes the simplicity of using the preassembled operators, although Definition (52) seems dense at first glance. Note that the stress localization depends on the temperature-dependent interpolation matrix φ (α) and on the coefficients of the secant interpolation.
The affine structure of the operator K (α) and of the right-hand side F (α) is closely related to the reduced basis method developed by Maday et al. [28] with ξ being the vector of boundary conditions and φ (α) being the sought after quantity of interest (see also, e.g., [32]). The use in the mechanical homogenization context and the reformulation of the mechanical problem using the general polarization strain (11) has not been exploited to the best of the authors' knowledge. Further, the use of distinct temperature intervals with local validity of the surrogate model resembles local reduced basis methods (see, e.g., [16]).

Greedy sampling
The accuracy of the proposed optimal interpolation method depends, as with any other interpolation method, on the width of the considered temperature interval and on the nonlinearity of the approximated quantities. Hence, an efficient sampling strategy is sought-after. We opt for a scheme that requires as few DNS evaluations as possible while allowing for an accurate prediction of the homogenized response.
An equidistant sampling of the temperature over the full range might not be the best approach. Hence, an enrichment indicator that is independent of the loading direction is sought. The relative nodal force error is evaluated seven times corresponding to ζ β ← e β (β = 1, . . . , 7), and the mean of these errors define an a posterior enrichment indicator e f that can be efficiently evaluated via where R •β refers to the βth column of the matrix R, • ∞ is the maximum norm and B ε is the standard FE strain operator. The normalization factor r ref is used to allow for a direct comparison with DNS residuals; it comes from the ground truth problem and is defined as corresponding to six orthonormal macroscopic strains ε γ = B γ (γ = 1, . . . , 6).

Numerical examples
Here, two materials that are used in laser melt injection are considered as RVE constituents: copper and fused tungsten carbide. Their temperature-dependent material parameters are approximated via a cubic polynomial with the form Polynomial coefficients corresponding to the utilized material parameters (visualized in Fig. 2) are given in Table 2. Most noteworthy is the decay of Young's modulus of copper toward 1300K where the contrast with respect to Young's modulus of fused tungsten carbide is 65.07 compared to 3.13 at room temperature.

Error measures
As error measures, we consider the average relative nodal force error Eq. (54) and the volume average of the point-wise error of the stress localization operator, in Eq. (18), evaluated via It measures the accuracy with respect to the reconstruction of the local stress field, which is deemed to be of major interest in potential applications. In addition, core outputs in thermomechanical homogenization are the effective quantities, namely the effective material stiffness C and the effective thermal strain ε θ . Therefore, corresponding error measures are defined next. Since C is symmetric and positive definite, its Cholesky decomposition is used to construct e C , i.e., the reference stiffness is written as We propose the two relative error measures: where I ∈ R 6×6 is the identity matrix.

Comparison between the proposed approaches
In this example, the accuracy of the proposed approaches is evaluated, and a comparison with a naive strain field interpolation is provided. In order to have a microstructure with an anisotropic response yet with a simple structure, a striped laminate structure is considered as the first example.
Virtual tests are run on 100 equidistant temperatures in the range [300, 1300] K. Consequently, a phase contrast of Young's modulus ranging from 3.13 at 300 K to 65.07 at 1300 K is observed which exemplifies the strong temperature dependence of the underlying homogenization problem. Empirical cumulative distribution curves (Fig. 4) use direct numerical simulation data at the two extremal temperature θ (1) = 300 K and θ (2) = 1300 K, i.e., no intermediate data are used in this example.
The first error e f corresponds to the evaluation of the relative nodal force error. Figure 4 illustrates the cumulative distribution function of e f . Ideally, the residual should be as close as possible to zero (meaning a jump of the plot from zero to one at x ≈ 0% is aspired). Then, the approximated solution will have a negligible error in the high fidelity model and, hence, accurate predictions of the underlying problem, i.e., Eq. (10). The response of approach (O 4 ) is basically indistinguishable from the reference solution, while the maximum error of 14.6% is observed for naive interpolation (N). Note that as expected and motivated earlier, approach (O 4 ) is the most reliable approach with 99% of the scenarios having a relative error smaller than 7 × 10 −13 %. (Such error is on the order of numerical roundoff errors.) At the same ratio of 99%, approaches (N, O 0 , O 1 , O 2 ) show relative errors up to 11.4%, 8.03%, 5.45% and 0.29%, respectively. Most notably, O 2 and O 4 are distinctly superior to the remaining approaches, immediately motivating their potential usefulness. The error e S is evaluated next. The error is visualized in Fig. 5 where it is seen that approaches (O 2 ,O 4 ) behave very well which indicates that approximating each of the material stiffness and the thermal strain separately is beneficial and leads to considerable improvement in approaches (O 2 ,O 4 ) over the results of approach (O 1 ). Furthermore, approach (O 4 ) guarantees the best approximation of the material stiffness and the thermal strain. It introduces relative stress localization errors (e S ) that do not exceed 5 × 10 −13 %, i.e., the method yields exact local stress predictions.

Illustration of the greedy sampling approach
Virtual tests are continued on another simple microstructure shown in Fig. 6. The temperature range stays the same as in the previous example [300, 1300] K. It is sampled via 100 equidistant points and the thermomechanical response of the RVE is evaluated on 100 uniformly random strain directions with a constant amplitude of one.
Now, e f is used to enrich the initial samples at 300 K and 1300 K in a hierarchical manner by successively adding further sampling temperatures one by one (each requiring 7 linear problems to be solved using DNS): The peak value of e f over all candidate temperatures is used to select the next sampling temperature. By construction, the error is numerically zero at the sampling points which is nicely seen in Fig. 7, which reports  the corresponding relative nodal force error e f . It is noted that adding two sampling temperatures (i.e., using data at four temperatures overall) reduced the relative error from more than 4.5% when using the min./max. temperatures only to less than 0.66%.

Comparison with an artificial neural network approach
A comparison with a black box artificial neural network (ANN) is presented in this example. The model uses the temperature θ as the sole input feature and predicts the effective quantities C and ε θ . Since machine learning is not the main topic of this work, interested readers are referred to, e.g., Alzub et al. [1]; Bishop and Nasrabadi [2] for further details about machine learning.
For the comparison between the proposed O 4 approach and ANN, three RVEs with randomly oriented, and shaped cuboids and spheroids are generated with various inclusion volume fractions of 20%, 40%, and 60%. These RVEs are visualized in Fig. 8 Remark 1 Note that while O 4 offers full-field reconstruction capabilities of all mechanical fields on the RVE scale, such reconstruction would be demanding using machine-learned models such as ANNs. Hence, the presented ANN model predicts only effective quantities, and the comparison with O 4 will be confined to C and ε θ .
Results corresponding to the first RVE on the left of Fig. 8, i.e., the one with 20% inclusions, are illustrated in Fig. 9. Figure 9a, b shows the superiority of O 4 that led to errors less than 0.02% and 0.2% for approximating Remark 2 The ANN model showed high approximation errors with a very low number of samples. Hence, the comparison with O 4 is confined to the best scenario with six samples.

Remark 3
The supposedly easy ML-centric approach-though known to potentially fail-does not even get sufficiently accurate predictions on the effective quantities, not mentioning the local fields at all. This is due to the limited amount of training data. It is worth noting that in the current application and the black box framework, more than 50 out of 100 data points were required for the training and validation of the ANN model to reach an approximation error smaller than 1%.

Staggered model for practical two-scale applications
The homogenized response consisting of C and ε θ is often sufficient for practical applications. As we have demonstrated, ANNs could be used for this purpose but they would require more input data. However, since our proposed scheme (O 4 ) allows the rapid precomputation of an abundance of data points with high accuracy, we can also rely on a staggered model. In such a model, trivial linear interpolation is used to obtain the homogenized response at arbitrary intermediate temperatures based on sufficiently resolved data.
This staggered model with space-saving tabular data has many advantages over the previously presented ANN approach, e.g., no model training or hyperparameter tuning is required, and accuracy can be increased arbitrarily by adding new samples using (O 4 ). An interface to this procedure is given in the GitHub repository [35] accompanying this article. With this approach, our presented scheme can be conveniently used to provide the homogenized response for practical two-scale applications. Moreover, this approach can also be used to generate the high amount of training data that other ML models (ANN, kernel regression, …) require for accurate and reliable predictions.

Conclusion
An efficient optimal field interpolation approach is proposed and investigated in this work. It is shown that such an approach has a very low numerical cost. Its cost is comparable to linear interpolation, but it is superior when it comes to accuracy. In many scenarios, the predicted quantities match their reference values. Optimal field interpolation is used to predict full thermomechanical fields of various RVEs with temperaturedependent materials. DNSs are required and evaluated at a few discrete temperatures. Then for any intermediate temperature, an energy optimal basis is derived from the sampled DNSs. Furthermore, a comparison with a black box machine-learned model is presented where the limitations of the ML model are illustrated in view of the limited availability of training data.
Applications of optimal field interpolation include but are not limited to: • Approximation of thermal homogenization problems that can be seen as a special case of the presented thermomechanical one, • efficient generation of training data for deep material networks [27] or other ML models, • and for sure an extension to plasticity will show the versatility of the proposed method.
A key advantage for the latter will be the access to local stress and strain fields offered by the proposed method.
When targeting the generation of tabular data (see also Sect. 6.5), one can start with a few calls to the high-fidelity finite element (FE) model in order to build the proposed ROM. Then, the ROM allows for rapid computations of the effective quantities over a very wide input range, i.e., for the generation of an abundance of high-fidelity data. The availability of many samples of such quality data is ideal; for instance, when it is fed into an ANN model, the ANN will benefit from the large amount of training data to ensure its accuracy. Hence, our ROM allows ML models to be trained efficiently and accurately. As illustrated in Sect. 6.5, the availability of large amounts of quality data can also eliminate the ANN modeling (with training, etc.) by replacing it with simplistic look-up table interpolation that can still yield arbitrary accuracy needs.

Supplementary information
The proposed interpolation algorithm is available from Sharba et al. [35]. All presented results may be reproduced using this repository. It is also possible to use the proposed approach on new RVEs given that their DNSs are stored in a HDF5® file that has a matching structure to the files we use Sharba and Fritzen [34]. Further details are available in the readme file in Sharba et al. [35].

Conflict of interest
The authors declared that they have no conflict of interest.

Appendix A: Mandel notation
Symmetric second-order tensors such as stress σ or strain ε are identified by their corresponding sixdimensional column vector, in Mandel notation, i.e., The previous expression is a result of utilizing the following symmetric orthonormal basis tensors of second order [7], where {e 1 , e 2 , e 3 } is an orthonormal basis of the three-dimensional Euclidean space E 3 . Fourth-order tensors with minor and major symmetries, e.g., C ↔ C ∈ R 6×6 are represented as where C •••• are the elements of C.