A numerical study of different projection-based model reduction techniques applied to computational homogenisation

Computing the macroscopic material response of a continuum body commonly involves the formulation of a phenomenological constitutive model. However, the response is mainly influenced by the heterogeneous microstructure. Computational homogenisation can be used to determine the constitutive behaviour on the macro-scale by solving a boundary value problem at the micro-scale for every so-called macroscopic material point within a nested solution scheme. Hence, this procedure requires the repeated solution of similar microscopic boundary value problems. To reduce the computational cost, model order reduction techniques can be applied. An important aspect thereby is the robustness of the obtained reduced model. Within this study reduced-order modelling (ROM) for the geometrically nonlinear case using hyperelastic materials is applied for the boundary value problem on the micro-scale. This involves the Proper Orthogonal Decomposition (POD) for the primary unknown and hyper-reduction methods for the arising nonlinearity. Therein three methods for hyper-reduction, differing in how the nonlinearity is approximated and the subsequent projection, are compared in terms of accuracy and robustness. Introducing interpolation or Gappy-POD based approximations may not preserve the symmetry of the system tangent, rendering the widely used Galerkin projection sub-optimal. Hence, a different projection related to a Gauss-Newton scheme (Gauss-Newton with Approximated Tensors- GNAT) is favoured to obtain an optimal projection and a robust reduced model.


Introduction
Phenomenological constitutive models are frequently used to compute the material response of a continuum body. However, the main influence of the macroscopic response is driven by the heterogeneous microstructure, whereas the direct modelling of the underlying microstructure is usually not feasible. A variety of analytical methods exists to account for the microstructure, e. g. Eshelby [12] or Mori-Tanaka-Method [27]. Since these methods are often limited, for instance to the linear regime or by the shape of the inhomogeneities that can be modelled, a different approach is necessary for the general nonlinear case and arbitrary shapes of the inhomogeneities. One possibility is given by computational homogenisation, which requires a nested solution scheme involving the computation of a boundary value problem (BVP) at the microscopic level, using a representative volume element (RVE) for every so-called material point.
Using the Finite Element Method to compute approximate solutions to the governing equations on both scales is often referred to as the FE 2 -method and does not rely on macroscopic constitutive models, but on the solution of underlying BVPs and a consistent transition between the two scales [16]. This approach is usually accompanied by high computational costs, due to the repeated solution of numerous BVPs at the micro-scale in a possibly high dimensional space. One approach to lower the computational cost is given by methods that rely on the fast Fourier transformation (FFT) [28,37]. These methods replace the assembly of the system tangent and the subsequent solution of the linear system on e.g. the micro-scale by an algorithm that utilizes FFT. This lowers the complexity of the solution process but does not necessarily reduce storage requirements [14].
Within this work, projection-based ROM is considered, which is characterised by taking advantage of the observation that the solutions of the aforementioned boundary value problems often lie in a lower dimensional subspace and different computational tasks in the offline and online stages. During the offline stage all necessary computations aiming to build a basis for the reduced-order model are carried out. This basis is used in the online stage to compute approximate solutions using a lower dimensional description. Several methods for model reduction exist, e.g. Proper Generalized Decomposition [8,10], Reduced-Basis methods [3,30,31], or approaches using the Proper Orthogonal Decomposition (POD). The latter has been widely used, e.g. in homogenisation, fluid mechanics and many other fields [6,17,18,23,26,41]. For a deeper discussion and a broader overview the reader is referred to [2,33] and references therein.
In the context of computational homogenisation the Reduced Model Multiscale Method (R3M) has been presented in [41], applying POD-based model order reduction to the BVP on the micro-scale for the case of hyperelasticity, directly projecting the governing equations onto the truncated POD basis. Due to the missing approximation of the arising nonlinear terms the computational savings were limited. In general, the approach of solely applying e.g. Galerkin projection for nonlinear problems may even lead to higher computational cost compared to the Finite Element model.
A further contribution in the context of geometrically linear homogenisation is given by [18], in which the stress field itself is approximated by a POD basis. Furthermore the authors showed in which circumstances an ill-posed system is obtained in the context of computational homogenisation and possibilities to avoid this problem. A similarly approach, i.e. not approximating the displacement field, is used in the nonuniform transformation field analysis (NTFA) [24], based on [11], and its extensions [14,15,25]. These methods apply a decomposition of the internal variables using reduced bases and derive suitable evolution equations for the reduced variables.
A popular approach in projection-based model reduction is the application of a Galerkin projection, as for e.g. the R3M, to solve the arising overdetermined system for the reduced coordinates. In the presence of nonlinearities these problems are often handled using an interpolation technique [1,7] or Gappy-POD reconstruction [5,13], both also often referred to as hyper-reduction techniques. Convergence difficulties have been reported for certain ROM configurations applying Galerkin projection coupled with hyper-reduction, e.g. [5,9,22,34]. In [9] the condition number of the reduced system tangent is suspected to lead to said convergence problems. The authors therefore propose a gappy data reconstruction with Galerkin projection to improve the condition number. While improving the robustness the number of diverging ROMs has only been reduced. In [5] an alternative to the Galerkin projection is proposed, which is related to the Gauss-Newton algorithm (GNAT). This projection lowers the constraints of the arising system tangent matrix and renders a robust model reduction technique.
The aim of this contribution is to apply projection-based ROM techniques in the context of computational homogenisation of hyperelastic materials including hyper-reduction techniques, which requires a robust reduced model for the present multi-query context. The focus thereby lies on the problem on the micro-scale, i. e. the quantities computed on the micro-scale which would be used on the macro-scale in a fully coupled problem formulation. The main objective is to compare three different model reduction approaches in terms of accuracy, robustness and optimality.
The remainder of this contribution is organised as follows: Sect. 2 discusses the fundamentals of first-order computational homogenisation and the governing equations. Section 3 describes the ROM techniques used within this study, including considerations regarding optimality of different projection techniques. Section 4 presents numerical examples, followed by some concluding remarks in Sect. 5.

Computational homogenisation in a nutshell
In solid mechanics, phenomenological constitutive models are frequently used to describe the material response of a body under a given load. However, the response is mainly driven by the heterogeneous microstructure. One possible approach to account for the heterogeneous material would be to directly model the substructure, e. g. within a FE model, which usually leads to computationally expensive models. A different approach is given by computational homogenisation. In this context the material at the macroscopic level is modelled without a prescribed constitutive model. Instead, the constitutive response is computed for every so-called material point at the micro-scale, taking into account the heterogeneities through prescribed constitutive behaviour of the constituents. To distinguish between the two scales, the superscripts (•) M and (•) m are used to denote quantities on the macro-and micro-scale respectively. Figure 1 illustrates the general concept of this method.
On both scales the balance of linear momentum represents the equation of interest, which reads on the micro-scale with F m representing the micro-scale deformation gradient and Ψ m the strain energy density. The material response on the macro-scale relies on quantities computed on the micro-scale. In order to compute these quantities, certain requirements have to be satisfied.

Hill-Mandel condition
A main ingredient of scale transition is the equality of the virtual work on the macro-and micro-scale, known as Hill-Mandel condition [19][20][21], where V 0 denotes the volume of the RVE in the reference configuration. Assuming a linear ansatz for the deformation on the micro-scale, withũ denoting the fluctuation field and X m , x m the coordinates in the reference and spatial configuration respectively, the deformation gradient is given by Considering the variation of the deformation gradient δF m = δF M + ∇ X δũ (2.6) and inserting Eq. (2.6) into Eq. (2.3) leads to It follows that the volume average of the microscopic Piola stress has to equal its macro-scale counterpart, i. e.
for the first term to vanish. Furthermore, the second term, containing the fluctuations, has to vanish in order to comply with the Hill-Mandel condition. Reformulating this term as a boundary integral renders admissible boundary conditions. Within this work the fluctuations are set to vanish on the Dirichlet boundary, which yields linear displacement boundary conditions: It should be noted that alternatively periodic boundary conditions could be applied.

Computation of the tangent
To provide information about the behaviour at a material point due to an increased load, the macroscopic tangent modulus has to be computed. In the present work the approach in [39] is used. Therein, a series of numerical test cases is performed in order to compute the fourth order Lagrangian elasticity tensor This leads, depending on the dimension d of the problem, to a total of d 2 linear perturbation test cases at the converged equilibrium state with the perturbations of the macroscopic deformation gradient with e i and E J denoting the basis vectors in the spatial and the reference configuration respectively. Hence, the computations are necessary to compute the constitutive constants.

Weak form and spatial discretisation
Within the scope of this contribution, an approximate solution to the micro-scale problem is computed using the FEM, which requires the weak formulation The discretisation of the above equations is carried out using the Bubnov-Galerkin Finite Element Method in a standard manner, i. e.
where the displacement field and the test function are approximated as a sum of shape functions (Lagrange polynomials) and nodal values, with n en denoting the number of element nodes. Using the definitions given in Eq. (2.14), the discretisation of Eq. (2.13) takes the form which may be written as a vector-valued equation Hence, f m represents a nonlinear function of the unknown nodal displacement values. In order to find an approximate solution an iterative Newton-Raphson scheme is used. This requires the linearisation of the nonlinear function which yields the tangent stiffness matrix Starting at an initial guess of the solution u m,k , an iterative update of the displacement is computed via leading to This iterative procedure is repeated until a prescribed convergence criterion is satisfied. Since the focus lies on the micro-scale problem, the superscripts (•) m will be omitted in the following in order to improve readability and only be used if it is required in the context. Furthermore, Eq. (2.17) will be solved for the unknown fluctuation fieldũ instead for the micro displacement field u.

Reduced-order modelling
As presented in the previous section, one possibility to compute approximate solutions to the microscopic problem defined by Eq. (2.1) involves the use of the Finite Element Method. Depending on the discretisation of the considered domain this may lead to large-scale systems. Especially in a multi-query context, as is the case in computational homogenisation, this leads to high computational costs. The solutions of such systems often lie in an affine subspace of lower dimension and therefore techniques to reduce the dimensionality of such problems are desired. Within the scope of this study a POD-based ROM approach is applied, including various hyper-reduction techniques.

Proper Orthogonal Decomposition
Consider discrete values a (•) j , e.g. displacement fluctuation values of the BVP of the micro-scale problem, solved using the Finite Element Method. These so-called snapshots are arranged into a matrix A with rank d ≤ min(n, n s ), n denoting the number of degrees of freedom of the FEM model and n s the number of snapshots. These solutions span a certain space denoted by V . The snapshot POD [36] is then used to filter out the dominant characteristics, allowing the computation of an orthonormal basis, that best suites a rank l d approximation of the snapshots in a least-squares sense. This task can be formulated as a constraint optimization problem. It can be shown that the solution of this optimization problem is given by the first l left singular vectors U (:, 1 : l) of A (ũ) = U · Σ · V T called the POD basis of rank l [2,40]. The basis vectors optimally represent the snapshots in a least-squares sense for the given rank l approximation. For the choice of a suitable l, it is useful to consider a truncation criterion ε in order to select the first l POD modes. For a system of rank d the criterion may be defined in terms of the singular values σ i as which gives information about the ability of the truncated basis to reproduce the snapshots. In the following, POD bases will be abbreviated by U r (•) , e.g. U r (ũ) ∈ R n×l for an ldimensional POD basis of the displacement fluctuation field u. In case of large snapshot sets, a nested POD as given in [4,32] may be used, which in essence partitions the snapshots into smaller sets, computes a lower rank approximation of each set and eventually computes a POD of the low rank approximations of the snapshot sets.

Projection approaches
Introducing the dimensionality reduction for the primary unknown via the POD renders an overdetermined system of equations and therefore suitable projection techniques are required. The widely used Galerkin projection and an alternative Petrov-Galerkin projection are briefly reviewed in this section. For a more detailed discussion the reader is referred to [5,35].

Galerkin projection
Employing the Galerkin projection the fluctuation fieldũ and the test function are approximated using the POD basis vectors, i. e.  (3.4) and the corresponding reduced tangent stiffness matrix Recalling the dimension of the POD basis matrix U r (ũ) ∈ R n×l , with n equal to the number of degrees of freedom and l the number of selected modes according to a chosen error criterion, the obvious benefit of this procedure is that now a system of equations for only l unknowns has to be solved. This decreases the computational cost especially for l n. Using the iterative Newton-Raphson scheme leads to the updated approximate solution As shown in [5,35] this projection, which produces Eq. (3.6), is optimal in the sense that it minimizes the error between the solutions of the reduced and the full order model in the K-norm: Thereby the tangent stiffness matrix has to be symmetric positive definite. It will be shown in Sect. 3.3 that neither of the hyper-reduction techniques discussed in this work guarantees the symmetry of the unreduced tangent stiffness matrix given in Eq. (3.20). Hence, the Galerkin projection combined with the hyper-reduction approaches as discussed within the work is not optimal in the sense of Eq. (3.8).

Petrov-Galerkin projection
As shown in Eq. (3.10) This approach renders and corresponds to the least-square problem requiring the tangent stiffness matrix solely to be regular [5,35]. While reducing the number of unknowns, the nonlinear terms still have to be evaluated at the full scale and projected onto the subspace at every iteration step, which clearly limits the computational savings. Hence, further reduction techniques have to be applied in order to significantly reduce the computational cost.

Hyper-reduction
As previously highlighted, the direct projection approach still depends on the full scale dimension n, due to the evaluation of the nonlinear terms. There exists a variety of approximation techniques for nonlinearities such as Empirical Interpolation Method (EIM) [1], its extension Discrete Empirical Interpolation Method (DEIM) [7] or the Gappy-POD [5,13], amongst others. It should be noted that, as shown in [18], employing hyper-reduction may lead to ill-posed systems, since the internal force vector, which is approximated using hyper-reduction, is zero at a converged state. Within our studies the bases for the subsequent hyper-reduction techniques are computed using snapshots of the internal force vector during the iterative solution process (the vector is non-zero). Hence, using a non-truncated basis for the approximated nonlinearity, one obtains the same internal force vector as that of the full order model, which justifies this approach. Within the present study the DEIM and the Gappy-POD in combination with the discussed projection approaches are compared and will therefore be shortly discussed.

Discrete empirical interpolation method
In essence, this method approximates a nonlinear function as where the parameter μ is introduced to denote the dependence of the fluctuation field on the macroscopic deformation gradient F M , used to compute the macroscopic displacement field. The parameter stems from a suitable parameter space μ ∈ D ⊂ R d , e.g. μ = F M 11 , F M 12 , F M 21 , F M 22 ∈ D ⊂ R 4 , for the two dimensional case. The direct projection approach requires the collection of snapshots a (ũ) i of the Finite Element approximated fluctuation field arranged into A (ũ) in order to compute the POD basis. Using DEIM, snapshots of the corresponding nonlinear function a (f) i , see Eq. (2.17), are collected during the iterative solution procedure in the offline phase and assembled into A (f) , where n s equals the number of considered snapshots. Performing the POD of A (f) renders the matrix U r (f) ∈ R n×k , representing a k-dimensional orthonormal basis, i.e. k modes are considered, for the space spanned by the snapshots of the nonlinear term. The coefficients of c in Eq. (3.13) are computed using k rows of f (ũ (μ)) where P denotes an extraction operator. This may be considered as a matrix composed of k vectors where i ρ i = [0, ..., 0, 1, 0, ..., 0] T denotes a vector in which the position of the only nonzero entry corresponds to the index ρ i [7]. Since the matrix P T · U r (f) is always regular [7] the coefficients of c can be uniquely determined. This leads together with Eq. (3.15) to the DEIM approximation The nonlinear term f now only needs to be evaluated at k entries specified by P. The corresponding DEIM indices ρ are determined using algorithm 1, proposed in [7], which computes the indices ρ based on the basis U r (f) . The reduced nonlinear term reads after Galerkin projection where the first term, of dimension l × k, represents a constant quantity and is thus computed during the offline phase. Online only k components, corresponding to the k DEIM indices, need to be computed. The tangent is obtained as the derivative of Eq. (3.18), The part of (3.19) which represents the tangent approximation, may not be symmetric as pointed out by [5,34]. Hence, applying a Galerkin projection is not optimal in the sense of Eq. (3.8). The same holds for the Gappy-POD in combination with a Galerkin projection. Consider therefore the following short example with the quantities and the sampling matrix Using these matrices to evaluate Eq. (3.20) produces which is not symmetric. This small example shows that it can not be guaranteed that the tangent approximation in Eq. (3.20) preserves symmetry.

Gappy-POD
Contrary to the interpolation in Eq. (3.17) the Gappy-POD [5,9,13] uses regression to approximate the nonlinear function. The approximation of the nonlinear term results in Here, k s indicates the number of sampling points with k s ≥ k, i.e. more sampling points than modes (keeping in mind that U r (f) ∈ R n×k ) and † denotes the pseudo-inverse. The tangent is computed analogously to Eq. (3.19). Similar to [5,9], algorithm 2 represents the point selection algorithm used in this work.

Gauss-Newton with approximated tensors (GNAT)
Instead of combining Galerkin projection and Gappy-POD, the GNAT solves the least-square problem in Eq. (3.12) using a Gappy-POD approximation of the nonlinear terms, which reads with the matrices T · U r (K) = I ∈ R k×k , while being independent of the dimension of the FEM model n. Here, the quantity U r (K) is introduced to account for the possibility of different snapshot selection strategies for the gappy approximation of the residual and the system tangent as presented in [5,6]. Within the scope of the present work snapshots of the residual obtained from the Finite Element model (including the iterative states during the Newton-Raphson solution procedure) are used. These serve as the input to build the reduced basis for both the residual and the tangent, i.e. U r (K) = U r (f) .  Based on this training set a few POD basis modes for the fluctuation field and the nonlinear function are depicted in Fig. 3. One may observe the influence of taking snapshots of f during the iterative Newton-Raphson procedure where the internal force vector is non-zero.

Robustness considerations
In this section the aforementioned different model reduction approaches, i.e. DEIM, Gappy-POD and GNAT are tested for various dimensionalities of U r (ũ) and U r (f) using the test set which is different to the set D train . Therefore each configuration is tested against 1296 testcases. In case that any of said test cases does not converge using the reduced model the configuration is highlighted with an "x" in the subsequent result plots. Otherwise the color indicates the relative error of the fluctuation field in the L 2 norm  Figure 4 shows the results for the DEIM approach as described in Sect. 3. One may observe that for various configurations the method rendered scenarios where at least one test case did not converge. In the multi-query context this clearly is an undesirable result as the model should return the quantities of interest for arbitrary input parameters. Considering the errors, they behave as expected and decrease for an increasing number of basis vector U r (ũ) . The influence of the dimension of U r (f) appears to be rather small after a certain threshold, see also Figs. 8 and 9.
In [9] it has been shown that using a Gappy-POD instead of pure interpolation benefits the condition number of the reduced tangent and should lead to a more robust model. This has also been observed within this work as shown in Fig. 5. While rendering a more robust model with respect to changes of the input parameters, as well as more accurate results, there were still configurations that lead to diverging test cases. This might be due to a possible loss of symmetry through the application of a hyper-reduction method and the subsequent non-optimal Galerkin projection. To illustrate the asymmetry ofK a comparison of the sparsity patterns for the given example from Fig. 2 is shown in Fig. 6.  Using GNAT to compute D test for various dimensions of the reduced bases with k s /k = 2. A "x" denotes the case were at least one test case within D test did not converge Employing therefore the GNAT, suited for non-symmetric tangent matrices, lead to the most robust model within the scope of this study as shown in Fig. 7. One can observe that no configuration lead to a diverging test case and the method appears to be very robust. Furthermore, equivalently for DEIM and Gappy-POD, increasing the dimension of U r (ũ) decreases the error, while a change of the dimensionality of U r (f) has only minor effect on the accuracy above a certain threshold. This gets confirmed considering Figs. 8 and 9 which depict the errors for the averaged stresses and the macroscopic tangent. Fig. 10 The reduced residual (r = −X · P T · f) from GNAT versus the relative error of the fluctuation field rel u with k s /k = 2; for every l the error is averaged over all test cases within D test Employing Gappy-POD with Galerkin projection rendered the most accurate results while the GNAT rendered the most robust model and showed no convergence difficulties. Though the Gappy-POD is more robust compared to the DEIM, there is no guarantee that a different D test would yield the same behaviour. It should be noted that the reduced model has been used for the perturbation procedure from Sect. 2.2 to obtain the macroscopic tangent components while the stresses were averaged over the full domain using the solution of the reduced model U r (ũ) ·û. A further benefit of GNAT is that it minimizes the global full order residual. Therefore using the reduced residual r = −X · P T · f from Eq. (3.23) may serve as an error indicator, considering the relation to the relative error as depicted in Fig. 10. The reduced residual as well as the relative error decrease with an increasing size of the reduced basis and match quite well up to a certain offset. This may be beneficial considering greedy selection algorithms to build the ROM, which rely on the ability to estimate or at least indicate the error of the reduced model approximation.
The relation between the errors rel u and the residual is further investigated for a given ROM and all test cases of the parameter set D test in Fig. 11. Based on the previous results the ROM dimensions are set to l = 45, k = 90 and k s /k = 2 using GNAT. It can be seen that the relative error has minor fluctuations but does not show considerable deviations from its mean. The same holds for the norm of the reduced residual from the GNAT model. Hence, the norm of the reduced residual appears to be a computational cheap and reliable error indicator which could be used in future studies in conjunction with an adaptive sampling algorithm, e.g. [32], instead of using D train to build the reduced model. Fig. 11 The norm of the reduced residual compared to the relative error of the fluctuation field for all test cases within D test using GNAT with l = 45, k = 90 and k s /k = 2

Local fields
Besides acting as a simple input-output model it is possible in computational homogenisation to investigate local fields. Considering a test case with

Conclusion
Within the scope of the present work reduced-order modelling techniques based on the Proper Orthogonal Decomposition and so-called hyper-reduction techniques have been applied in the context of computational homogenisation of hyperelastic materials. The focus has thereby been the accuracy and robustness of the reduced model, comparing different hyper-reduction and projection approaches. It has been shown that introducing an additional approximation of the nonlinear terms via an empirical interpolation or Gappy-POD may not preserve the symmetry of the system tangent. This leads to a non-optimal Galerkin projection. As shown in the numerical examples this can cause convergence problems. This is clearly an undesirable property in the multi-query context as given in computational homogenisation. The Gauss-Newton like approach (GNAT), relying on a Petrov-Galerkin projection suited for non-symmetric tan-gent matrices, rendered the most robust model and showed no convergence problems while minimizing the global residual. Future studies should aim towards more effective methods to obtain error estimates for the outputs of interest and averaging procedures. Furthermore, alternative approximations of the system tangent, e.g. [5,29,38], could be investigated. Yet some of these methods require additional high dimensional snapshots of the sparse system tangent which becomes computational infeasible rapidly. For instance considering only the non-zero elements of the tangent from the presented examples in Sect. 4 renders an additional snapshot of the size 40.836 for the 2.312 unknowns in every iteration step.
Regarding the snapshot selection adaptive strategies, e.g. [32], could be employed to adaptively select the positions in parameter space for which the full model should be evaluated in order to build the ROM.