A reduced order model for geometrically parameterized two-scale simulations of elasto-plastic microstructures under large deformations

In recent years, there has been a growing interest in understanding complex microstructures and their effect on macroscopic properties. In general, it is difficult to derive an effective constitutive law for such microstructures with reasonable accuracy and meaningful parameters. One numerical approach to bridge the scales is computational homogenization, in which a microscopic problem is solved at every macroscopic point, essentially replacing the effective constitutive model. Such approaches are, however, computationally expensive and typically infeasible in multi-query contexts such as optimization and material design. To render these analyses tractable, surrogate models that can accurately approximate and accelerate the microscopic problem over a large design space of shapes, material and loading parameters are required. In this work, we develop a reduced order model based on Proper Orthogonal Decomposition (POD), Empirical Cubature Method (ECM) and a geometrical transformation method with the following key features: (i) large shape variations of the microstructure are captured, (ii) only relatively small amounts of training data are necessary, and (iii) highly non-linear history-dependent behaviors are treated. The proposed framework is tested and examined in two numerical examples, involving two scales and large geometrical variations. In both cases, high speed-ups and accuracies are achieved while observing good extrapolation behavior.


Introduction
Driven by advances in additive manufacturing and tailorable effective properties of metamaterials, there has been a growing interest in understanding structure-property relationships of complex microstructures.These microstructures can typically be described by a few shape parameters, leading to distinct types of effective behavior.To investigate such structure-property relations and to find the optimal shape for a given application, simulations are often considered.These simulations are in general computationally expensive or even intractable for direct numerical simulation, especially for large-scale engineering applications, since considerably fine meshes are required to capture the complex microstructural geometry.By employing multi-scale methods based on computational homogenization [1,2] or domain decomposition methods [3], such large-scale problems can be separated into many smaller subproblems, thus rendering them amenable for efficient numerical simulation.Domain Decomposition (DD) methods are particularly useful when the micro-and macroscale are of comparable size, i.e., in the absence of a clear scale separation.They can be categorized in overlapping and non-overlapping DD methods.Regarding the latter, the domain is divided into subdomains and coupled at the interfaces.One notable method is the so-called FETI-DP [4], where the problem is solved in the corner points of each subdomain and in the Lagrange multipliers that enforce the interface continuities.To increase the computational efficiency, model order reduction methods were combined with DD.In the Reduced Basis Element (RBE) method [5], each subdomain is accelerated with a reduced basis and the interfaces are coupled weakly in a non-conforming manner with Lagrange multipliers.In [6], the Static Condensation Reduced Basis Element (SCRBE) was introduced where the internal degrees of freedom of each subdomain are represented by a reduced basis and condensed out, resulting in a conforming approximation space on the interfaces (also referred to as ports in this context).Constructing optimal local approximation spaces for these ports in two-component systems was discussed in [7], and finding them by local solutions of the Partial Differential Equation (PDE) with random boundary conditions was proposed in [8].In the context of solid mechanics, recent applications of such methods include, for example, [9,10].For a more comprehensive overview of concepts in localized model order reduction, the interested reader is referred to [11].
If scale separation is assumed, i.e., when the length scale of the typical microstructural features is much smaller than that of the macrostructure, first-order computational homogenization can be employed.Here, the behavior of the microstructure dictates the (average) constitutive behavior of an effective macrostructural continuum model.By defining a Representative Volume Element (RVE) which models the fine-scale geometry of the microstructure in full detail, a coarse-grained representation of the macrostructure with a much coarser discretization can be assumed at the macroscale.At every macroscopic integration point, the macroscopic strain is used to specify a microscopic boundary value problem which, after solution, returns the effective stress and stiffness.Since a PDE needs to be solved at every macroscopic Gauss integration point, this methodology is still computationally expensive, and efficient ways for its solution are needed.
Several approaches to tackle this problem have been reported in the literature.One class of methods are based on data-driven reduced order models.For these methods, the microscopic problem is solved several times to generate training data, and the data is subsequently used to learn an effective constitutive model.After training, the microscopic solver is not required anymore and replaced by the learned constitutive model.Several methods were proposed for elastic material models, see, e.g., [12][13][14][15], and also for historydependent microstructures [16,17].Even though highly efficient and accurate reduced order models can potentially be obtained, the extension of the methods [12][13][14][15] for inelastic materials is challenging, and the methods for history-dependent microstructures [16,17] often require vast amounts of data.To learn an elasto-plastic material model, Mozaffar et al. [16] generated up to 15000 deformation paths each of 100 load steps.9000 deformation paths of up to 2000 steps each were generated by Wu et al. [17].
Another class of methods attempts to accelerate the existing microscopic solver.For instance, if the Fast Fourier Transform (FFT) [18,19] is used to simulate the microstructure, its solution can be accelerated by the (nonuniform) Transformation Field Analysis (see, e.g., [20,21]), or Self-consistent Clustering Analysis [22,23].One disadvantage of FFT is that geometrical parameterizations of the RVE cannot be directly treated and, hence, sensitivities for material optimization cannot be directly computed.If the microscopic problem is solved via the Finite Element (FE) method, the resulting multi-scale formulation is referred to as FE 2 [24][25][26].By directly solving the microscopic PDE with FE, material or shape parameterizations can be considered in a straightforward manner, making the approach more suitable for inverse problems and optimization.To speed up the microstructural simulation, Proper Orthogonal Decomposition (POD) [27,28] can be utilized to find a reduced set of basis functions; the method then computes the Galerkin projection of the solution onto the space spanned by the snapshots.Although POD generally requires many full-order solves for training, it typically works well for all input parameters.In the context of first-order homogenization, POD was first applied in Yvonnet et al. [29] for a hyper-elastic RVE, and later explored in Radermacher et al. [30] for an elasto-plastic RVE under small strains.However, due to the nonlinearities of the microscopic problem, the speed-ups were limited since the global force vector and stiffness matrix must be assembled by full integration in every microscopic Newton iteration.To address this issue, a further reduction called hyperreduction is required, which aims at finding an efficient way of assembling microstructural force and stiffness quantities.Notable hyperreduction methods are Empirical Interpolation Method (EIM) (see, e.g., [31]), a variant of EIM called Discrete Empirical Interpolation Method (DEIM) (see, e.g., [32]), energy-based mesh sampling and weighting [33], reduced integration domain [34], empirical quadrature procedure [35], and Empirical Cubature Method (ECM) [36].EIM and DEIM interpolate the non-linear integrand of the global force vector such that the integrals can be pre-computed.In [37,38], DEIM was used successfully to accelerate the solution of the microscopic PDE.However, these works only discussed the solution of the microscopic PDE and did not derive the effective stress and stiffness quantities required for the macroscopic problem.A possible disadvantage of EIM and DEIM is that they lead to nonsymmetric tangent matrices, which might result in convergence issues, observed in, for instance, [38,39].The rest of the above-mentioned hyperreduction methods aim at approximating the integrals by finding a subset of integration points with corresponding positive weights among the set of all integration points used in the formulation of the microstructural PDE.This has the advantage that the stiffness matrix is always symmetric and at least positive semi-definite (in practice usually positive definite unless instabilities occur), ensuring a good convergence of the microscopic problem.Example applications of hyperreduction methods have been successfully employed in two-scale simulations in [40], where an elasto-plastic composite RVE under large deformations was considered.In [41], a damage model for a composite RVE under small deformations was shown.While both works obtained accurate results and successfully accelerated the forward simulations of a two-scale problem, such formulations were limited to fixed microstructures only, i.e., did not account for possible parameterizations.In order to allow for optimization of microstructures, the surrogate model needs to be extended to a wide range of different design parameters (including geometrical as well as material).
This work aims to address this gap, by developing a hyper-reduced surrogate model for geometrically parameterized microstructures, to enable (shape) sensitivity analysis and optimization of materials.Furthermore, we intend to provide a detailed analysis of the reduced RVE problem for arbitrary loading paths and geometries and elucidate possible issues due to reduction.Our main contributions are: 1. development of a hyper-reduced POD model for a family of geometrically parameterized microstructures, by employing a geometrical transformation method [42] and by extending the ECM algorithm to geometrical parameters, 2. consistent derivation of the effective stress and stiffness of the hyper-reduced model, 3. an empirical analysis of the accuracy of the surrogate model for elasto-plastic RVEs under large deformations for different geometries and loading conditions, 4. a quantitative comparison of a two-scale example with continuous change in microstructural heterogeneities.
The remainder of this paper is organized as follows.In Section 2, the microscopic problem arising in first-order computational homogenization is briefly summarized, together with the computation of the effective quantities.Section 3 covers in-depth development of the reduced order model with particular focus on the empirical cubature method for geometrically parameterized microstructures, and includes a detailed derivation of the effective stress and stiffness.In Section 4, the proposed method is examined and tested in detail, first for a single RVE and then also for a full two-scale problem.Finally, a summary on the findings and concluding remarks are given in Section 5.
In this work, the following notational convention is adopted.Italic bold symbols are used for coordinates X and vectorial or tensorial fields, such as the displacement field u or stress field P .Upright bold symbols are used for algebraic vectors and matrices, such as the global stiffness matrix K or the coefficients of the discretized displacement field u.A field quantity u for given parameters µ is denoted as u(X; µ).Given second-order tensors A and B, fourth-order tensor C, and vector v, the following operations are used: where the Einstein summation convention is implied.

Formulation Of The Microscopic Problem
In multiscale schemes based on first-order homogenization, the macroscopic problem is governed by the standard linear momentum balance, but the macroscopic constitutive model (relating strains to stresses and stiffness) is replaced by a microscopic PDE (again governed by the standard linear momentum balance) which is defined on an RVE.By prescribing the macroscopic deformation gradient on the microscale, the PDE can Figure 1: Two-scale problem based on first-order homogenization.At every macroscopic point, a microscopic simulation is defined through deformation gradient F and shape parameters µ, and solved to obtain an effective stress P and stiffness Ā.For different macroscopic points, different parameterized microstructures can be considered through µ.As an example of a family of geometrically parameterized microstructures, a parent domain with a circular inclusion Ω p (center), can be mapped onto parameterized domains Ω µ 1 (left) and Ω µ 2 (right) with mapping Φµ 1 and Φµ 2 .
be solved and an effective stress and stiffness returned to the macroscopic solver, see Fig. 1.For applications such as microstructure optimization, it is reasonable to additionally introduce a parameterization of the RVE in order to compute the sensitivities with respect to design variables.The microscopic boundary value problem is formulated below on a parameterized domain, as is usually the case in shape optimization.For brevity, the dependence on the macroscopic coordinates is omitted and a fixed macroscopic material point is assumed unless otherwise specified.

Boundary Value Problem
Consider a family of domains Ω µ ⊂ R d with space dimension d = 2, 3, parameterized by geometrical parameters µ, and spanned by a position vector X µ ∈ Ω µ .In Fig. 1, an example parent domain with a circular inclusion Ω p is geometrically parameterized and mapped to two distinct parameterized domains with elliptical inclusions, Ω µ1 and Ω µ2 .The volume and the topology of the domain |Ω µ | are assumed to remain fixed for all parameters (the outer boundaries of the RVE domain are fixed while the shape of the interior geometry can change).With the assumption of scale separation between macro-and microscale, the microscopic displacement field on the parameterized domain u(X µ ) can be written as the summation of a mean field ū(X µ ) and a fluctuation field w(X µ ), i.e., u(X µ ) = ū(X µ ) + w(X µ ).The mean field is fully specified through ū(X µ ) := ( F − I)X µ , where F is the macroscopic deformation gradient tensor and I is the identity tensor.The total deformation gradient tensor F is defined as The governing microscopic PDE is given as where Div(•) is the divergence operator with respect to X µ and P denotes the second-order first Piola-Kirchhoff (1PK) stress tensor.No constitutive model is specified at this point, although we assume that the stress P is a non-linear function of the deformation gradient F (or its history).The weak form of the problem is then: given the macroscopic deformation gradient F , find the fluctuation field where the integral bounds depend on the parameters µ, δw denotes a test function, and H 1 (Ω µ ) is a Hilbert space with square integrable functions and square integrable derivatives.The inner product in V is defined as From Eq. (3), it is apparent that the macroscopic deformation gradient F represents the external loading, while the fluctuation displacement field w balances the system.To simplify the problem in Eq. ( 3) and remove the parameter dependence of the integral bounds, a parent domain Ω p is defined.To this end, we assume that there exists a parameter-dependent diffeomorphism Using integration by substitution, the problem of Eq. ( 3) can be restated as follows: given the macroscopic deformation gradient F , find with the transformation gradient The superscript p is used to denote quantities pertinent to the parent domain, e.g., w(X µ ) = (w • Φ µ )(X p ) = w p (X p ).To iteratively solve the non-linear problem in Eq. ( 5), a linearization using the Gateaux derivative around the current state w p in direction ∆w p ∈ V p is required and can be written as, where A := ∂P ∂F is the fourth-order stiffness tensor.Once the transformation map Φ µ is known, Eq. ( 3) can be solved on the parent domain using Eqs.( 5) and (6).Further details on how to find these transformations for a range of geometrical parameters are provided in Section 3.4.By following a standard Galerkin finite element discretization for w p ≈ w p h ∈ V p h ⊂ V p , with dim V p h = N , the number of degrees of freedom of the dicretization, the internal force vector f ∈ R N and global stiffness matrix K ∈ R N ×N can be derived from Eqs. ( 5) and ( 6), resulting in the following non-linear system of equations where w ∈ R N is the column vector of unknown coefficients of the discretized fluctuation field.This problem can be solved with the Newton method, i.e., where m is the Newton iteration number and Eq. ( 8) is repeated until ||f (w m )|| 2 ≤ ε newton with ε newton a user-defined tolerance.For more details on the finite element method and discretization of weak forms, we refer to [43].

Effective Quantities
For conciseness of notation, the following abbreviations are introduced to denote quantities after the solution w * p has been obtained: Upon obtaining solution w * p from Eq. ( 5), the effective stress is computed as and the effective stiffness (in index notation) as where I mnkl := δ mk δ nl is the fourth-order identity tensor.To determine , Eq. ( 5) is differentiated with respect to F .For one particular component Fkl (where the indices k and l are assumed to be temporarily fixed), the differentiation yields where a new auxiliary vector field q kl := ∂w * p ∂ Fkl ∈ V p has been defined (reflecting the sensitivity of the microfluctuation field with respect to the change of the applied macroscopic loading), and E kl ∈ R d×d is a second order tensor with all entries zero, except for the kl-th entry which is 1.The linear tangent problem of Eq. ( 13) is then solved for all combinations k, l = 1, ..., d to obtain q kl for each component of F .
Although not utilized in this work, the sensitivities of the effective stress P with respect to the geometrical parameters µ, which are required for applications such as shape optimization, can be computed with the geometrically parameterized formulation of the RVE as follows (in index notation), The integrand is complicated due to the derivatives of F −1 µ and | det F µ |, but in principle these derivatives can be computed for a given geometrical mapping Φ µ .If the effective stress P can be assumed to vary smoothly with the parameters µ, which may be a reasonable assumption for smoothly varying shapes using, for instance, splines, finite differences can be used to approximate these sensitivities.

Surrogate Modelling
Since the microscopic problem has to be solved at every macroscopic quadrature point, the solution of the microscopic PDE must be efficient.Solving it directly using FE is in general too computationally expensive, and, hence, the microscopic solver must be accelerated.In this section, a surrogate model for the geometrically parameterized microscopic PDE is developed by employing the Reduced Basis Method (RBM) [27] to reduce the number of degrees of freedom and the Empirical Cubature Method (ECM) [36] to reduce the number of quadrature points.The key idea is to construct the surrogate model on the parent domain Ω p , adapt it to each geometry Ω µ , and then solve the reduced problem.

Reduced Basis Method
For complex problems and geometries, typically a fine mesh is required for FE, leading to a highdimensional solution space V p h for the fluctuation displacement field w p with dim V p h = N .The idea of the RBM is to approximate the field with global parameter-independent basis functions and parameterdependent coefficients, i.e., where N is the number of basis functions, ideally much smaller than the dimension of the FE space, i.e., N ≪ N .The basis functions, {ϕ n } N n=1 , span a subset of V p h and can be obtained by applying proper orthogonal decomposition (POD) on a set of pre-computed full solutions for different parameter values.Additionally, they are orthonormal with respect to V p , i.e., where δ mn denotes the Kronecker delta.By utilizing the POD space for both the trial and test space and inserting w p from Eq. ( 15) into Eqs.( 5) and ( 6), the components for the reduced internal force vector f POD ∈ R N and reduced global stiffness matrix K POD ∈ R N ×N can be derived as where a = [a 1 , . . ., a N ] T is the column vector of unknown coefficients to be solved for, and i, j = 1, . . ., N span over all basis functions.Analogously to Eqs. ( 7) and ( 8), the resulting non-linear system of equations can be solved using Newton method:

Empirical Cubature Method
Even though the solution field and linear system of equations have been reduced to dimension N ≪ N , computing the components of the force vector in Eq. ( 17) and global stiffness matrix in Eq. ( 18) still requires integrating over the RVE.For the full integration, a numerical quadrature rule (usually based on Gauss quadrature) with integration points and corresponding weights {( Xq , ŵq )} Q q=1 , where Q is the total number of integration points, is employed, i.e., for i = 1, . . ., N .For a fine mesh, Q is very large and thus evaluating Eq. ( 21) leads to high computational costs.To address this issue, we employ the Empirical Cubature Method (ECM), which was proposed in Hernández et al. [36] for a fixed geometry, and extend it to parameterized geometries.
The idea of ECM is to find a subset of points set of all integration points with corresponding weights {w q } Q q=1 that approximates Eq. ( 21) up to a user-defined error ε.To find such a subset that approximates Eq. ( 21) well for all admissible geometrical parameters µ, Eq. ( 21) is first rewritten as where the weighted stress W is defined.To remove the parameter dependence of the integrand in Eq. ( 22), the weighted stress is approximated by another reduced basis, i.e., where {B l } L l=1 is a set of L basis functions obtained using POD, which are orthonormal with respect to L 2 (Ω), i.e., Inserting Eq. ( 23) into Eq.( 22) and rearranging yields, Since Eq. ( 25) should be accurate for any choice of coefficients α l ( F , µ), all the N • L terms in Eq. ( 25) that approximate the integral have to be approximated as accurately as possible.Hence, the goal becomes to find a subset Q(≪ Q) of integration points with corresponding weights {(X q , w q )} Q q=1 that approximates Eq. ( 25) well, i.e., These Q points and corresponding weights are found using a greedy algorithm, the details of which can be found in [36] and are omitted here.The algorithm is terminated when the mean squared error of all N • L terms is less than a user-defined tolerance ε.
Compared to the original algorithm for a fixed geometry, as proposed in [36], the only differences are that the weighted stress W is employed instead of the stress P and that the parent domain Ω p is considered instead of a fixed domain Ω.With the ECM integration rule, the hyper-reduced force vector and global stiffness matrix are computed as Remark 3.1.The computational costs of the ECM greedy algorithm as proposed in [36] increase drastically with the number of selected integration points, since for every selected point a non-negative least squares problem needs to be solved.As pointed out in [44], rank-one updates can be used with the least squares solver for better efficiency, and such a refined version of the ECM algorithm was presented in [10].For the numerical examples considered in this work, the original ECM algorithm in [36] was sufficiently fast and we did not use the algorithmically improved version.

Effective Quantities
Once the new set of integration points and weights is found, the integrands of Eqs. ( 27) and ( 28) only need to be evaluated at the points {X q } Q q=1 during the solution of the reduced problem.This also means that the stress and stiffness field are available at these points only.To compute the effective quantities, the most straightforward method is to use the integration rule obtained by ECM, i.e., P = |Ω p | −1 Since the stress field P * p is known at all integration points {X q } Q q=1 , the effective stress can be directly evaluated.The method yields very accurate results in the examples considered below in Section 4.However, it should be noted that there is currently no guarantee that the integration rule found by ECM will generally be accurate for the computation of the effective stress.In general, the effective stress can have two sources of error as compared to the full solution: one comes from the solution of the reduced system and one from an inaccurate integration of the obtained stress field.The ECM integration points are selected such that the first error is minimized, but this also indirectly affects the second one to decrease, although not as quickly.This can be observed in the results of the first numerical example presented in Section 4.1.To ensure an accurate integration of the effective stress, it could be included into the ECM algorithm as a criterion.
As discussed in Section 2.2, derivatives ∂w * p ∂ F are needed to find the effective stiffness Ā, see Eq. ( 12).For each component of F , the linear tangent problem of Eq. ( 13) needs to be solved.By employing the trial space of the fluctuation field for the auxiliary function q kl , i.e., and the integration rule found by ECM, the following linear system of equations results: where q = [q 1 , . . ., q N ] T is the column vector of unknowns to be solved for and Note that the matrix K * p ∈ R N ×N is exactly the same as the global stiffness matrix K of Eq. ( 28) evaluated at the solution w * p .After solving the tangent problems, the effective stiffness Ā can be computed (in index notation) as where

Auxiliary Problem For Geometrical Transformation
Thus far, the geometrical transformation Φ µ : Ω p → Ω µ has been assumed to be known and has not been discussed in more detail.However, such transformations are in general not known analytically and have to be found numerically by using, for example, radial basis functions, see, e.g., [45,46], or mesh-based methods, see, e.g., [42,47].For each of those methods, an auxiliary problem arises which needs to be solved.In order to rapidly solve the surrogate model for a wide range of different geometries, it must therefore also be ensured that the auxiliary problem can be solved rapidly.In this work, the method in [42] is employed, in which the auxiliary problem is formulated as a linear elasticity problem by defining Φ µ (X p ) = X p + d(X p ), with d the transformation displacement obtained from Div C aux : 1 2 In the above equation, C aux is the fourth-order elasticity tensor, fully specified by the Young's modulus E aux and Poisson's ratio ν aux .The boundary conditions for this PDE are problem-dependent and are specified by the geometrical parameters µ.For the RVE problem, the outer boundaries are fixed (d = 0), while d is prescribed on parts of the interior that are parameterized by µ.In [42] the effect of the choice of E aux and ν aux was studied and it was demonstrated empirically that the choice only has a minor effect on the final approximation quality.Hence, in all numerical examples considered in this work, a Young's modulus of E aux = 1 and Poisson's ratio ν aux = 0.25 is assumed.The auxiliary problem can then be significantly accelerated with the RBM in combination with a (D)EIM [31,32], resulting in where Â ∈ R Np×Np is the reduced system matrix, d ∈ R Np is the reduced transformation displacement, b(µ) ∈ R Np is the reduced forcing vector and N p is the number of geometrical parameters.Since N p is usually small, Eq. ( 37) can be rapidly solved.From d, the transformation gradient F µ , its inverse F −1 µ , and its determinant det F µ can be computed.Moreover, expressions for the derivative of the inverse and determinant of the transformation gradient F µ can be derived, which are needed for computing the sensitivities with respect to µ, see Eq. ( 14).For more information on the auxiliary problem and its reduction, the reader is referred to [42].

Summary
For convenience, the offline-online decomposition for constructing and solving the surrogate model is summarized in Algorithm 1.
Algorithm 1: Offline-online decomposition of the proposed PODECM framework with microstructures parameterized with external loading F and geometrical features µ.
Offline Stage: 1: Define a parent domain Ω p and its finite element discretization.
2: Generate parameter samples { F i , µ i } Ns i=1 from a random distribution.3: For each different set of geometrical parameters µ i , solve the auxiliary problem in Eq. ( 36) to obtain the transformation map Φ µ i .
4: Compute F −1 µ i and det F µ i for each parameter sample µ i , then run full simulations (Eqs.( 5) and ( 6)) for F i and collect fluctuation displacement and weighted stress snapshots.5: Compute POD for the fluctuation displacement and weighted stress, cf.Eqs. ( 15) and ( 23).6: Run ECM algorithm and find integration points and weights, cf.Eq. ( 26).7: Assemble the reduced system matrix and forcing vector for the auxiliary problem in Eq. ( 37) by applying POD and DEIM.Details are provided in [42].Online Stage: 1: Given a new parameter set ( F * , µ * ), solve reduced auxiliary problem Eq. ( 37) and compute F −1 µ * and det F µ * .2: Solve reduced problem for F * with Eqs. ( 27) and ( 28).3: Compute effective stress using Eq. ( 29).4: Solve the linear problem Eq. ( 31) for each component of F * .5: Compute components of the effective stiffness with Eq. (34).

Example Problems
The proposed framework, referred to as PODECM, is first tested on a non-linear composite microstructure under various loading conditions and analyzed in depth regarding its capabilities and accuracy.The RVE consists of an elasto-plastic matrix with stiff inclusions of variable size and is considered under nonmonotonic loading.The surrogate model is analyzed in terms of the number of basis functions of the fluctuation displacement field N , number of basis functions of the weighted stress L and the ECM integration error tolerance ε.Subsequently, a two-scale problem involving a porous microstructure under non-monotonic loading conditions and varying porosities is studied to illustrate the accuracy and speed-up of PODECM in a two-scale setting.
All experiments are defined in two dimensions under plane strain conditions.The RVEs are assumed to be of size [0, 1] 2 and all quantities are assumed to be normalized and hence dimensionless.Since the macroscopic deformation gradient F can always be decomposed into a rotation R and a symmetric stretch tensor Ū with a polar decomposition, i.e., F = R Ū , it is sufficient to generate training data for the stretch tensor Ū , having only 3 independent components (6 in 3D).
To measure the quality of the approximation, the following error measures to compare the full FE simulations against PODECM solutions are defined: 1. Error of effective stress where P PODECM ( Ū k ) and P FE ( Ū k ) denote the effective stress obtained with PODECM and FE for Ū k , || • || F denotes the Frobenius norm, K is the total number of loading steps and Ū k is the applied external load at load step k.

Error of fluctuation field
where w PODECM and w FE denote the fluctuation displacement field obtained with PODECM and FE, and .f., Eq. ( 4).Recall that the integral in Eq. ( 4) is defined over the parameterized domain Ω µ .

Elasto-Plastic Composite RVE With Random Stiff Inclusions 4.1.1. Problem Description
The considered RVE in this example consists of two phases, an elasto-plastic matrix and stiff elastic inclusions.The geometry of the parent domain is shown in Fig. 2a, where the volume fraction of the inclusions is 23.4%.For the geometrical parameterization, one geometrical parameter µ = {ζ} that scales the size of the inclusions uniformly (and is proportional to the volume fraction of the inclusions) is introduced, see Figs. 2b and 2c showing two example domains for distinct values of ζ.The simulation mesh is depicted in Fig. 2d, where six-noded quadratic triangular elements are used in conjunction with three quadrature points per element.In total, the mesh has 62194 degrees of freedom, 15450 triangular elements and 46350 quadrature points.For the constitutive model of both matrix and inclusion the small-strain J 2 -plasticity model with linear isotropic hardening is chosen and extended to large strains with the method presented in Cuitino and Ortiz [48].The details of the plasticity model are provided in Appendix A. For the matrix, the following material parameters are selected: a Young's modulus E = 10, Poisson's ratio ν = 0.3, yield stress σ y0 = 0.2 and hardening constant H = 5.For the inclusions, E = 100 and ν = 0.3 are selected, corresponding to a stiffness contrast ratio of 10 between both components.Since no plastic deformation is assumed for the inclusions, their yield stress is set to a large value such that yielding never occurs. Three

Results
In total N train = 20 samples are generated from a Sobol sequence to train PODECM whereas 100 testing samples are generated from a uniform distribution to test PODECM.Each sample consists of 40 snapshots for each load step.
The accuracy and speed-up of PODECM depends on the number of basis functions used for the fluctuation displacement field N and the number of quadrature points Q.While N is typically chosen directly, Q depends on the choice of the number of basis functions used for the weighted stress L and the ECM integration error ε.
To study the influence of L on the resulting number of quadrature points Q and mean errors in effective stress and fluctuation field on the testing dataset, several combinations of N and L for a fixed ε = 10 −2 are tested, with resulting errors shown in the top row of Fig. 4. The projection error (for N basis functions and using full integration) is shown as well.It can be clearly seen that the number of quadrature points Q increases drastically with increasing N and L, as more information needs to be integrated accurately.In fact, a roughly linear relationship Q ∝ N L can be recognized.For the mean errors, a higher L leads to better results on average, although we observe that errors fluctuate significantly, and for some values of N a worse approximation is obtained with a higher L. This occurs since the ECM algorithm is a greedy algorithm, meaning that it does not necessarily find an optimal set of integration points.When more basis functions are included into the algorithm, a completely different set of points may be found that finally leads to a worse approximation.It can furthermore be observed that the gap between the projection error and the PODECM solution grows larger for increasing N .This is because the basis functions typically become more oscillatory and difficult to approximate with higher N , see, e.g., [13,49], and thus require significantly more quadrature points for a good approximation.It is interesting that the gap for the errors in the fluctuation field are smaller than the ones in the effective stress, i.e., the difference between the computed PODECM and the projection errors are much higher for the effective stress as compared to the fluctuation field.This happens because the ECM integration points and weights are primarily selected to integrate the weak form accurately, and using them to compute the effective stress introduces an additional approximation error, cf.Section 3.3.
Several combinations of N and ε for a fixed L = 15 are next tested to study the influence of ε on the number of quadrature points and approximation errors.Obtained results are shown in the bottom row of Fig. 4. Similarly to the previous analysis, a lower ε leads to more quadrature points Q and a lower mean error in the effective stress and fluctuation field on average, as the integrals are approximated more accurately.Interestingly, lowering the tolerance from 0.01 to 0.001 does not significantly improve the approximation quality, even though substantially more quadrature points are included, meaning that the errors can be attributed to the higher modes of the weighted stress (the additional quadrature points barely contain any information).Therefore, choosing a tolerance smaller than ε = 0.01 leads to no improvement.From Fig. 4 we further observe that the errors of the fluctuation field are considerably higher (order of magnitude) than the errors of the effective stresses.This results from the fact that the POD basis functions aim to minimize the H 1 (Ω p ) error, and thus approximate the field accurately on average rather than locally, suggesting a favorable approximation for averaged quantities such as the effective stress.In [50], the authors showed a similar result for FFT: the effective stresses converge with an order of h, while the local fields converge with h 1/2 , where h is the used voxel edge-length.
To test the data efficiency of PODECM, the reduced model is trained for different numbers of training data N train = {1, 3, 5, 20} with L = 20 and ε = 0.01.The number of integration points and corresponding errors in the effective stress and fluctuations are shown in Fig. 5.It can clearly be seen that the number of integration points barely changes for different N train , and that the errors converged already for N train = 3.No noticeable improvements with N train = 5 and 20 can be observed, showing that PODECM is very data efficient.Even with N train = 1, the errors in the effective stress are below 5%.
To conclude, the more basis functions N and L are used and the lower the integration error ε is chosen, the more accurate the final result is.However, at the same time the surrogate model grows in size and the speed-up decreases.A user must thus make a compromise between accuracy and cost.In our experience, the speed-up correlates nearly linearly with Q, i.e., if the number of quadrature points is reduced by a factor of 100, this results in a speed-up of roughly 100.In contrast, the number of basis functions N only plays a minor role for the speed-up.For this example, use of N = 20, L = 20 and ε = 0.01 lead to a reduction in the number of degrees of freedom from 62194 to 20 and quadrature points from 46350 to 212, suggesting a speed-up on the order of roughly 200.
Remark 4.1.Intuitively, one might consider not selecting L as an independent parameter, but as a (nonlinear) function of N with L ≥ N , since the weighted stress acts as a non-linear function on the fluctuation field, and thus is expected to require more basis functions.However, for L = N , one obtains a roughly quadratic relationship for the number of quadrature points Q and N with Q ∝ LN = N 2 .If a higher number of N is necessary to have a sufficiently large solution space, Q quickly becomes very large and speed-ups of PODECM become very small.By treating L as an independent parameter and allowing L < N , the resulting number of integration points can be controlled and decreased.Furthermore, we observed in the numerical tests that, if a maximum number of integration points is specified, we often obtained better results for N > L rather than N ≤ L, as long as L is large enough.

Problem Description
In the second example, the macroscopic structure, depicted in Fig. 6 together with the employed simulation mesh, is compressed under an external loading T (x).Here, we assume H = 1, W = 2, with T the magnitude of the applied load.The simulation mesh has 1322 degrees of freedom, 200 8-noded quadrilateral elements, and in total 800 quadrature points.The structure is assumed to have a porous microstructure, modelled by the parameterized RVE, shown in Fig. 7, and the same material model as the matrix material in the previous example.Such microstructures (with circular holes, i.e., a = b) have been considered in several works, see, e.g., [51][52][53][54], due to their auxetic behavior under compression, i.e., negative Poisson's ratio.During compression, the center part of the material starts to rotate, thus pulling the material from the sides inwards.In our work, we define two independent parameters, namely the volume fraction of the voids, v void := 4πab, and the ratio of the semi-major axis b to semi-minor axis a of each hole, κ := b/a.The semi-minor axis a and semi-major axis b depend on v void and κ as    Depending on the values of the parameters, the resulting effective properties change significantly.To illustrate this, linear analyses of this RVE for different parameters have been carried out, similarly to [55], where a small compression in the y-direction with ∆u y = 0.001 has been applied, while allowing the RVE to contract freely in the x-direction.With the resulting displacement in the x-direction, ∆u x , the Poisson's ratio in the initial state can be estimated as Similarly, the initial Young's modulus is estimated as where Pyy is the yy-component of the effective stress.
The parent mesh is selected with κ = 1.25 and v void = 0.45.The external loading is applied in K = 50 load steps.In the first 25 load steps, the applied load T increases linearly from 0 to 0.2.In the next 25 steps T is decreased linearly from 0.2 to 0.

Results
Several two-scale simulations with different PODECM surrogate models are run and compared to the full reference FE 2 solution.To compare the accuracy of the surrogate solutions, the compliance C := Γ T (x)u y (x)dx, where Γ denotes the top horizontal edge of the macrostructure and u y its vertical displacement, is computed at every load step.The compliance is an important quantity, often employed in optimization problems.Subsequently, the relative error in compliance ϵ C and the relative error averaged over all load steps εC are defined as where the subscript k denotes the k-th load step and C FE 2 is the compliance computed with the full solution.
The tested PODECM RVE models are generated for different numbers of training samples N train and number of basis functions N .The training data is sampled from Ūxx ∈ [0.85, 1], Ūyy ∈ [0.85, 1], Ūxy ∈ [−0.15, 0.15], v void ∈ [0.4,0.5] and κ ∈ [1.01, 1.5] with a Sobol sequence.Each sample consists of 40 load steps, where the sampled macroscopic stretch tensor is applied to the RVE with a piecewise linear amplitude function that is linearly increased from 0 to 1 for the first 20 load steps and then linearly decreased from 1 to 0 for the last 20 steps.For the ECM algorithm, the number of basis functions for the weighted stress L and the integration error ε are assumed as fixed with L = 20 and ε = 0.01 for all models.The exact settings for each surrogate model are summarized in Table 1, alongside the averaged relative error in compliance εC , and the run time and online speed-up in comparison to the full FE 2 solution.For comparison, the total number of degrees of freedom and quadrature points of the full FE model of the RVE are also provided.As can be seen, for all surrogate models the number of quadrature points are reduced by a factor of up to 100 which results in high speed-ups up to 100 times, while errors are below 5% for all models.By increasing N from 10 to 50 for rom 1 to rom 5, the error decreases from 4.74% to 1.54%, whereas the speed-up reduces from 92 to 26.Including more training samples with the same number of basis functions N = 50 (from N train = 20 for rom 5 to N train = 50 for rom 6) improves the error from 1.54% to 0.39% while the speed-up remains roughly the same.This means that by increasing the sample size, the first 50 basis functions contain more general information that result in a better approximation.When more than 50 samples are used for rom 7 and rom 8, the error remains roughly the same, implying that 50 samples are sufficient for this problem and including more training data does not improve the results.In Fig. 9 the force-displacement curve is shown for the FE 2 and a few selected surrogate solutions.The displacement is defined as the vertical displacement at the mid point of the top edge (which is also the maximal displacement).It can be observed that all surrogate models underpredict the displacement, indicating that the surrogate models overpredict the stiffness of the macrostructure.
The relative error in compliance is plotted over the load steps k in Fig. 10.For all models, the error slowly increases over k.The reason for this behavior is that all training samples are generated for simple loading cases, where the macroscopic stretch tensor is linearly varied in only one direction throughout the entire simulation.In the macroscopic simulation, the macroscopic stretch tensor for one integration point generally does not evolve along one direction, but changes its direction continuously, leading to highly complicated deformation paths and histories that are not included in the training data.To tackle this problem, random loading paths during training could be used, as performed in, e.g., [16,17], to generate a more general surrogate model.Solely increasing the number of samples from 20 to 50 (rom 5 to rom 6) also decreases the observed errors to less than 1% for all load steps.This shows that PODECM generalizes well to loading paths that are not part of the training data.
Regarding the offline run times (for one core of an Intel Platinum 8260): computing each training sample takes approximately 2-3 minutes.With the training data available, computing the POD for each model takes up to around one minute.Finally, selecting the ECM points takes up to roughly two minutes for the PODECM models with N = 50.6a.The structure starts deforming plastically for T > 0.07 and a residual displacement of roughly ũ = 0.04 remains after unloading.All surrogate models achieve accurate results, although they generally predict a slightly stiffer response than the full FE 2 solution.

Conclusions
In this work, we developed a reduced order model, termed PODECM, by combining the proper orthogonal decomposition (POD), the empirical cubature method (ECM) and a geometrical transformation method.This method is designed to accelerate microscopic computations within two-scale simulations, parameterized by material and geometrical parameters, hence enabling PODECM to be used for two-scale optimization problems.We showed how to compute the effective stress and the corresponding consistent effective stiffness for this model, as required by the macroscopic solver, and how to obtain the effective sensitivities with respect to geometrical parameters.The framework was first tested on a single-scale problem involving an RVE of a composite microstructure that consisted of a soft elasto-plastic matrix with stiff inclusions of variable size, controlled by a single geometrical parameter.With PODECM, the number of degrees of freedom and integration points was reduced to a fraction of the full FE model while maintaining a high accuracy in effective stress.The performance of PODECM was further evaluated for a two-scale simulation, in which a porous microstructure, characterized by two geometrical parameters, was considered.Both geometrical parameters were varied throughout the macrostructure, and depending on their values, the effective Poisson's ratio changed from positive to negative.For this example, different PODECM models were constructed with good accuracies (of errors less than 1%), while achieving speed-ups up to 100.Even though highly accurate and fast solutions were obtained with the proposed method, several open questions prevail.For instance, optimality and accuracy of the ECM integration rule cannot be ensured and thus needs to be further analyzed and understood.
As the underlying microscopic PDE is still being solved, only a small amount of training data is required to construct a good approximation of the full model.History-dependent material behavior does not need to be specially treated, making the proposed framework very general and viable for two-scale shape optimization problems.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.1) begin with a lower error which slowly grows with k.By increasing the sample size of the training data (comparing rom 5 and rom 6), the prediction becomes more accurate for the same number of basis functions; the errors for all load steps for rom 6 are below 1%.

Figure 2 :
Figure 2: Parent with two parameterized domains and simulation mesh.(a) The parent domain consists of a matrix material (blue) with 23 random elliptical inclusions (orange).The problem has one geometrical parameter ζ that scales all ellipses uniformly (ζ = 1 for parent domain).(b) A parameterized domain for ζ = 1.2 and (c) for ζ = 0.5.(d) The considered mesh consists of six-noded triangular elements and contains in total 62194 degrees of freedom, 15450 triangular elements and 46350 quadrature points.

Figure 4 :
Figure 4: The left column shows the number of quadrature points Q obtained from ECM for different choices of number of basis functions of fluctuation field N , number of basis functions of weighted stress L, and ECM integration error ε.The middle and right columns show the average errors of the effective stress and the fluctuation field when tested on the testing data for different choices of N , L and ε.The top row assumes a fixed ε = 0.01, while the bottom row assumes a fixed L = 15.

1 εwFigure 5 :
Figure 5: Errors of PODECM for different numbers of training samples.The number of quadrature points does not change.For N train = 3, the errors in effective stress and fluctuation displacement is already converged, showing that PODECM is very data efficient.

Figure 6 :
Figure 6: Geometry (a) and mesh (b) of the considered macroscopic structure.The body is fixed on the bottom and an external compression force T is applied on the top.The mesh consists of 1322 degrees of freedom, 200 8-noded quadrilateral elements and 4 quadrature points per element.

Figure 7 :
Figure 7: Geometry (a) and mesh (b) of the porous RVE.The elliptical holes are all characterized by the same semi-minor axis a and semi-major axis b, and are parameterized by the volume fraction of the pores v void and the κ = b/a.The employed simulation mesh has 21042 degrees of freedom, 4964 6-noded triangular elements and in total 14892 quadrature points.The parent domain corresponds to κ = 1.25 and v void = 0.45.
For parameter ranges v void ∈ [0.4,0.5] and κ ∈ [1.01, 1.5], the estimated Poisson's ratio and Young's modulus are plotted in Fig.8.It can be observed that removing material (by increasing v void while keeping κ fixed) or increasing κ while keeping v void fixed both Effective Young's modulus E eff

Figure 8 :
Figure 8: Initial effective Poisson's ratio (a) and Young's modulus (b) of RVE for different values of v void and κ.

Figure 9 :
Figure9: Force-displacement curve of two-scale simulations.The displacement ũ is defined as the vertical displacement of the mid point on the top edge, cf.Fig.6a.The structure starts deforming plastically for T > 0.07 and a residual displacement of roughly ũ = 0.04 remains after unloading.All surrogate models achieve accurate results, although they generally predict a slightly stiffer response than the full FE 2 solution.

Figure 10 :
Figure10: Relative error ϵ C in compliance over load step k.All surrogate models (specified in Table1) begin with a lower error which slowly grows with k.By increasing the sample size of the training data (comparing rom 5 and rom 6), the prediction becomes more accurate for the same number of basis functions; the errors for all load steps for rom 6 are below 1%.

Table 1 :
Summary of results for full FE2and different PODECM surrogate models.Reduced order models are generated for different numbers of training samples N train and basis functions of the displacement N .The number of quadrature points Q follow from the ECM algorithm with a fixed number of basis functions of the weighted stress field L = 20 and an integration error of ε = 0.01.All reduced order models achieve errors less than 5% with speed-ups up to 100 times.By generating more training data and maintaining the same N (rom 5 and rom 6), better results are achieved.However, using more than 50 training samples (rom 7 and rom 8) leads to no further improvements, suggesting that 50 samples are sufficient for this problem.
a All computations are executed using 20 cores of an Intel Platinum 8260.