1 Introduction
Modern material design industry demand highperforming simulation tools capable of capturing the nonlinear effects of materials observed at the structural scale, although governed by the physics occurring at small length scales. During the last years, a considerable effort has been put into the development and improvement of multiscale computational homogenization techniques for analyzing heterogeneous composite materials Matouš et al. (2017).
Typically, for material modeling involving two wellseparated length scales, e.g., structural and mesoscopic scales, special attention has been paid to hierarchical techniques that show a high capability to capture the complex mechanical interaction effects between scales (Michel et al. (1999), Miehe et al. (1999), Terada and Kikuchi (2001), Blanco et al. (2016)). A pioneering paper of this kind of techniques is the fe2 approach by Feyel et al. Feyel and Chaboche (2000), which can tackle the nonlinear behavior of complex microstructures by solving surrogate problems on a rve of the target material Gitman et al. (2007); Nguyen et al. (2010).
However, fe2 approaches are still computationally unaffordable for their use in software for practical applications. This is the main reason why they have not yet been widely transferred to the industrial sector. As an example, an analysis involving a standard fe2 multiscale simulations may typically take from months to years to complete in typical computing clusters.
A general concept to defeat the barrier imposed by the overwhelming cost of fe2 techniques, without giving up the displayed accuracy to capture the interaction effects between scales, consists of separating the full computational burden of a multiscale material simulation in two stages. First, an offline training stage is carried out with the hf model of a rve, where a data set of hf solutions is sampled and a surrogate or reduced model of the homogenized material response is built. Finally, this reduced model, which demands a low computational cost, is used in the online structural virtual testing. Several methodologies that follow this concept are very briefly described:

Artificial neural networks can be used to construct a cheap parameterized surrogate model. Results reported in
Rocha et al. (2020) clearly show that the computational cost of these strategies during the online stage can be lower than the one demanded by alternative reduced models. However, surrogate models derived in this way ignore the physical basis of the problem, as the model parameters are nonphysical coefficients. Additionally, these strategies are unable to simulate material problems governed by input parameters unforeseen in the sampled load trajectories; e.g., for reproducing nonmonotonous loading conditions that were not encountered in the sampling. See Ghavamian and Simone (2019), where a technique to overcome this drawback in historydependent materials is proposed. 
Employment of wavelet functions (Deslauriers and Dubuc (1989), van Tuijl et al. (2019)
) to define different resolutions for the interpolation of basis functions. This strategy can be seen as an onthefly optimal integration technique for reduced models, at the expense of a low speedup—only one order of magnitude—with respect to a full rom.

The pgd method Chinesta et al. (2013) is a powerful reduction technique that has been used in a number of engineering applications. Within the pgd approach, a multidimensional solution is computed offline, where model parameters and boundary conditions are considered extra coordinates of the problem. The multidimensional solution is approximated by the sum of modes, composed by the product of the dimensional coordinates (including spacetime, material properties and boundary conditions). This product consists of a lightweight offline computational catalog to be employed for rapid simulations of industrial interest. Typical procedures adopting this approach have been addressed by Ladeveze and coworkers Ladevèze et al. (2010); Chinesta et al. (2011), Doblare and coworkers El Halabi et al. (2016), among others.

Hyperreduced finite element methods, resulting from a combination of a rom technique (which in turn is based on pod) with a procedure that reduces the computational complexity related to the nonlinear terms of the variational formulation. The pod technique provides an empirical lowdimensional basis of spatial interpolation functions that replaces the original fe interpolation functions of the hf model. Furthermore, several techniques to decrease computational complexity have emerged in last years. The objective is to select only a few Gauss integration points—the cubature points—where the nonlinear terms are computed. The parameters of these cubature points are associated with a physical term of the problem. Some examples of hyperreduction techniques are the discrete empirical interpolation Chaturantabut and Sorensen (2010), the energy conserving sampling and weighting method Zahr et al. (2017), the use of reduced integration domain Ryckelynck (2009), and the use of pod modes to approximate the integrand of the nonaffine terms, plus the selection of a reduced number of sampling points via gappy data reconstruction Hernández et al. (2014). A related technique to the present approach is the empirical cubature method reported in Hernandez et al. (2017), and improved in Hernández (2020).
Following the approach explained in point iv, the hprfe2 technique—described by the authors in previous contributions Oliver et al. (2017), Caicedo et al. (2019)—is evaluated in this paper. The principal characteristic of this model lies in its potential use as a software engine for the material design industry to overcome the tyranny of scales. It has been reported in LloberasValls et al. (2019) that this hprfe2 technique is able to provide a speedup of up to four orders of magnitude for microstructural discretizations of several millions of integration points, assuming relative errors of around with respect to hf solutions. These speedups imply that a multiscale simulations that originally could take up to several years, can now be performed in a few hours. This clearly indicates that hprfe2 technology fulfills the requirements of the material design industry for an affordable multiscale analysis tool.
It is emphasized that the hprfe2 formulation proposed in this paper differs from the Empirical Cubature Method proposed in Hernandez et al. (2017) in the sense that, in order to construct a reduced integration scheme, what is approximated is the fundamental (variational) principle, i.e., the free energy minimization, and not its derivatives (stresses and internal forces) as it will be described in the formulation section. The minimum of the variational principle (energy minimization) is, in any case, well approximated and since is a scalar field, the reduction procedure yields a lighter reduced basis. Consequently, the performance of the approach is considerably improved with respect to the approach from Hernandez et al. (2017)
, which is based on integration of vector entities. This is, among others, one of the specific features of this method.
The objective of this work is to further assess the potentiality of this technique in terms of a tradeoff between computational speedup, consistency, accuracy, and preservation of the physical basis of the model. The main novelty of the current contribution consists in the application of a HPRFE methodology to an industrial simulation case (fibrereinforced composite laminates). This entails a number of upgrades to the strategy with respect to former contributions: 3D analysis, tailormade constitutive modeling, full reconstruction of stress and internal variables fields, assessment of the real speedups and computation times for realistic simulation scenarios.
In the following Sections 2 and 3, a brief summary of the hprfe2 technique is presented (a detailed description can be found in Oliver et al. (2017), Caicedo et al. (2019), LloberasValls et al. (2019)). The accuracy and consistency of the methodology is assessed in Section 4, by comparing the solutions of two composites virtually tested with the hprfe2 and the hf techniques.
2 Multiscale problem formulation
Let a hierarchical multiscale material model such as the one schematized in Figure 1 be assumed. The coupon represents the macroscale structure, where effective mechanical variables are considered, i.e., displacement , strain , and stress . A rve of the microstructure is used to compute the homogenized constitutive response of the material at each point of the structure. Computational homogenization consists in transferring onto the rve, solving a microscale bvp, and upscale the effective stress evaluated as the volumetric average of the microstress
. The same operation also computes the corresponding effective constitutive tangent tensor
.The first point to be defined in the formulation is the admissible kinematics at the microscale. The displacement field at the microscale , results in the addition of three terms:
(1) 
where the last term is the displacement fluctuation field . Accordingly, the microstrain in the rve is given by the addition of two terms:
(2) 
where is assumed uniform in the microscale volume , while the microstrain fluctuation is the symmetric gradient of the displacement fluctuation field^{1}^{1}1In the remaining of this paper, the infinitesimal strain theory will be considered., satisfying:
(3) 
In accordance with this kinematics, a bvp at the rve is next formulated in terms of microstrain fluctuations Oliver et al. (2017).
MICROSCALE PROBLEM FORMULATION: Given the macroscale strains , and the spaces of kinematically compatible strain fluctuations and admissible strain fluctuations
(4) 
find such that
(5) 
(6) 
where the space of microstrain tensorial functions is defined such that its elements fulfil the infinitesimal strain compatibility conditions given by
(7) 
with being the dimension of the Euclidean space.
Equation (5) is the variational form of the equilibrium equation in terms of the microstress which, in turn, implicitly considers the constitutive relation as a function of the internal variable^{2}^{2}2In this work, a single scalar internal variable damage is considered. microscopic damage and microstrains. Conventional damage constitutive models are assumed at the microscale. Therefore, the damage rate equation is explicitly defined by (6).
The homogenized stress is the volumetric average of the microstresses which are the solution of the bvp (4)(6):
(8) 
and the homogenized constitutive tensor is:
(9) 
where are the constitutive tensor of the composite phases, is the fourth order identity tensor and is the localization tensor ().
Since the microscale bvp is formulated in terms of the strain field, the displacement does not play any role in the rve problem, and the history of the strain components are the only parameters determining the homogenized response^{3}^{3}3This is also a feature of the proposed formulation, that allows optimal specific treatment of zones of the rve with different constitutive equations..
By adopting (2) and a kinematically admissible microdisplacement fluctuation space satisfying the last identity in (3), the micromechanical bvp can be rewritten in a more conventional displacement formulation. As will be shown in the next section, both formulations, in displacements and strains, are used to develop the hprfe2 technique. The displacementbased formulation is used in the hf fe technique for the model sampling stage, and the strainbased formulation is employed within the reduced and hyperreduced models. The advantage of taking this approach is commented in the following section.
3 hprfe2
The hprfe2 model uses a lowdimensional space of functions for approximating the microfluctuation strain field , as described in Section 3.1. An identical approach is taken for the microstrain variations, i.e., a Galerkin formulation. Additionally, a reduced numerical cubature rule, the roec, is introduced in Section 3.2 to compute the integral balance in (5). Section 3.3 describes the procedure for recovering the strains, displacements and damage fields of the rve, i.e., projecting fields from the cubature points onto the original fe mesh.
3.1 Dimensional reduction of the microstrain fluctuation field
The lowdimensional space representing is built as the span of an orthogonal basis of spatial functions with global support: (where and is the dimension of the strain and stress tensors in Voigt notation for 3D problems), as follows:
(10) 
where each element of the basis is a microstrain fluctuation mode, and the vector of pseudotime dependent coefficients () represents the amplitude of these modes. In the last identity of (10), the matrix , with collects, in columns, the microstrain modes of the basis . Notice that, in order to preserve a simplified notation, an identical symbol identifies the microstrain fluctuation field in both, hf and lowdimensional approaches.
The basis , which is composed of basis vectors (or modes), is computed with a pod technique applied to a set of microstrain fluctuation fields, in turn obtained as solutions of the microcell problem (hf model formulated in displacements) during an offline sampling process.
Finally, the variational formulation (5)(6) is projected onto the space spanned by the basis . Therefore, the number of equations of the corresponding discrete nonlinear systems reduces to equations, whose solution determines the vector .
Remarkably, by construction, each element belongs to the vectorial space (7) inheriting the boundary conditions imposed to the hf model in the sampling stage. Therefore, any function spanned by the basis is an admissible microfluctuation strain.
3.2 roec
The integral term in the global balance equation (5) is computed with the roec rule. This rule is derived as follows (find additional details in Oliver et al. (2017)): the dimensionality of the free energy function space computed in the sampling stage is reduced using a similar approach to that adopted for the microstrain fluctuation field (10). Thus, assuming that the lowdimensional free energy is spanned by a basis of elements,
(11) 
The number of cubature points of the roec rule , their spatial position and the corresponding weights are selected with a similar criterion to that reported in Hernandez et al. (2017). This criterion is based on adopting an exact integration^{4}^{4}4Exact integration understood as the selection of the optimal set (among the Gauss integration rule adopted in the hf model) and its corresponding weights, in order to exactly integrate the selected energy modes . of the free energy basis , plus the condition that the reduced integration of the unitfunction in gives the volume of the microcell . These conditions are expressed as follows:
(12)  
(13) 
where the symbol refers to the Gauss quadrature rule in the hf model. This criterion provides conditions. Therefore, the number of cubature points which can be obtained is .
Remark: A remarkable feature of the roec technique hinges on the determination of the reduced integration rule in terms of the free energy field (fundamental variational principle) and not on its derivatives (stresses and internal forces), as proposed in many other formulations. As a result, the problem goal, i.e., find the minimum of the variational principle (free energy), is the actual key of the solution of the problem, unlike in alternative methods based on finding null values of the functional derivatives (internal forces)^{5}^{5}5I.e., the EulerLagrange equations of the variational principle.. Consequently, the formulation results simpler and very accurate, leading to the high performance reduction technique hprfe2.
3.3 Reconstruction of fields in the rve
After solving the hprfe2 problem, given the approximate microstrain field and the internal variables of the damage model ( described in A) in the cubature points, the microfields can be recovered and projected onto the original hf mesh as described below.
3.3.1 Reconstruction of the displacement fluctuation field
The displacement fluctuation fields can be recovered at the original nodes of the hf mesh by using the following procedure. Let the conventional space of the fe interpolation functions for displacement associated to the hf mesh be defined, but excluding the displacements of rigid body modes. Then, the displacement reconstruction problem is stated as follows:
FIND: such that:
(14) 
Introducing the fe approach of and the variations of microdisplacement fluctuations , and considering that is the global vector of nodal parameters, interpolating in the original fe mesh, then (14) can be rewritten as follows:
(15) 
where is the straindisplacement matrix of the conventional hf fe method (). Finally:
(16) 
Notice that matrix from (16) has to be computed only once for the reconstruction process.
3.3.2 Reconstruction of the internal variable field of the damage model
Similarly, the internal variable of the damage model (summarized in A) can also be recovered. Firstly, a lowdimension spatial approach for is defined. In the offline sampling of the original microcell, snapshots of this field are gathered and a reduced basis of elements is determined by means of a pod technique, similarly to the procedure defined in (10). Then, the reduced field can be expressed as:
(17) 
The vector (whose components are the internal variable values at the reduced cubature points and at a given pseudotime ) is obtained from the hprfe2 solution. Secondly, is projected onto the vector (being the total number of Gauss integration points), whose components are the internal variables at the Gauss quadrature points of the hf mesh. The following problem is defined:
FIND: such that:
(18) 
Applying a similar procedure to that used for the microdisplacement fluctuation recovery in Section 3.3.1, coefficients can be computed as follows:
(19) 
Finally, using equation (17), we recover the vector , whose th component (at the spatial position of the th Gauss point ) satisfies
(20) 
where . Since the damage variable is recovered by projecting a vector from a lower dimensional approximation space onto the higher dimensional space of finite elements, in this operation there may arise values lesser than 0 or greater than 1. Thus, it is necessary to disregard those values assigning to the recovered damage variable the value 0 or 1, respectively, to each case.
3.3.3 Reconstruction of the microstress tensor field
Given the internal variable vector and the microstrain at each Gauss quadrature point of the hf fe mesh, we compute the microstress tensor at the same points by appealing to the damage constitutive equation.
Remark: Most of the involved computations in (16) and (19) can be performed during the offline model sampling stage, typically, the evaluation of matrices and . The matrix vector products and , as well as the evaluation of the microstress at the Gauss integration points of the hf fe mesh can be computed during a postprocessing stage and after solving the hprfe2 problem. These additional evaluations represent a low computational cost.
4 hprfe2 model assessment
The hprfe2 technique is assessed by means of the multiscale simulation using two microstructures of engineering interest, representing glass fiberreinforced composites with epoxy matrix. Particular attention is paid to the tradeoff between attained accuracy, computational speedup, and the capacity of the hprfe2 model for simulating nonsampled trajectories.
The offline strategy of the sampling stage of both models is explained in Section 4.2, which details the number and directions of the loading scenarios needed to attain the reduced basis employed in different approximations of the hf solution. The error of the reduced order approximation with respect to the hf solution is calculated as explained in Section 4.3. Model completeness, i.e., the capacity for modeling material responses which have not been specifically considered during the sampling stage, is detailed in Section 4.4. Computational performance for microcells of increasing size and complexity is shown in Section 4.5. Material customization, i.e., the capacity of the model to reproduce the response of materials whose properties are different to those used during the offline model sampling, is evaluated in Section 4.6.
Examples of multiscale virtual tests are shown in Section 4.7, where a test coupon is simulated under varying configurations of the microstructure (obtained by rotations of the same microcell and customization of the material parameters), and a number of combinations of representation modes and integration points. Recovered damage and stress fields of the microcell in sampled elements of the coupon is also shown.
4.1 Microscopic models

Model A represents a periodic composite constituted by a generic crossply laminate containing several plies of aligned longitudinal fibers ( diameter, volume fraction of ) at and angles, embedded in an epoxy matrix. It represents a material that can manifest the modeled failure mechanisms (fibermatrix pullout, fibermatrix decohesion, and interply delamination), rather than a realworld material example. The rve of this composite, as well as the material properties of its constituent phases, are defined in Table 1. Fibers are modeled elastic, while matrix, interply, and fibermatrix cohesive interphase are assumed inelastic, modeled by a damage constitutive model with bilinear hardening, as described in A.

Model B represents one ply of an industrial multilayer composite. The microstructure consists of a random distribution of unidirectional fibers ( diameter, volume fraction of ) in an epoxy matrix. Fibers are assumed elastic in the strain range considered. Matrix is modeled by the same damage model referred previously. Contrarily to model A, the fibermatrix interaction effects ignore the fiber pullout and decohesion mechanisms. The material parameters of the constituent phases of this composite are reported in Table 2.
The objective pursued with the analysis of both models is not necessarily to perform their precise assessment and validation against experimental results. Rather, it is the demonstration of the capabilities of the hprfe2 technique for capturing the composite deformation modes involving discontinuous strain fields, particularly those arising in model A. Accurately capturing these deformation modes results in a challenging target for the present approach. To assess the hprfe2 capability, the solutions obtained with this methodology are compared to those obtained with the hf model, which are considered as the reference solutions in all tested cases.
Matrix  Fiber  Cohesive  Interply  
Young modulus []  
Poisson ratio  
Elastic threshold []  –  
Infinity elastic threshold []  –  
Hardening parameters ,  ,  –  ,  , 
Matrix  Fiber  
Young modulus []  
Poisson ratio  
Elastic threshold []  –  
Infinity elastic threshold []  –  
Hardening parameter ,  ,  – 
4.2 Sampling
The sampling stage is performed by solving a set of trajectories, i.e., problems in the rve with the hf technique, subjected to a macroscopic strain as the evolving action.
By considering that the six components of a symmetric macrostrain tensor can be related to a point lying on the space, we choose
points uniformly distributed in the unitary hypersphere in
Kunc and Fritzen (2019), to define the same number of corresponding macrostrains. Macrostrains are imposed to the rve multiplied by a time factor , monotonically increasing from 0 to for model A, and for model B. The coefficient is chosen such that most of the trajectories reach the inelastic regime, which assures a sufficient number of inelastic snapshots during the sampling stage. The solutions of the full loading trajectory set constitutes the rve sampling procedure. Although a damage model is employed to capture the material nonlinearities, only a moderate strain localization is expected since softening is not included in the constitutive response. Hence, the adopted sampling range suffices to correctly reproduce the nonlinear regime and it is not expected that severe localization takes place way beyond the one captured during the sampling stage. An extension of the method to account for softening would require a regularization strategy such as the one analyzed inOliver et al. (2017).The sampling data set is composed of snapshots (one for each timestep) per each trajectory, from where the pod bases are obtained. There are no special requirements for the selection of the snapshots, provided that the selected sampling trajectory represents the full parametric strain space. The approach taken here consists in gathering as many inelastic snapshots as possible, in a large trajectory diversity, and to leave the selection of optimal snapshots to the SVD stage. A snapshot consists of strain, energy and damage variables in every Gauss integration point (ip) of the hf mesh, at a given time step. These snapshots are required for building the reduced set of integration (cubature) points and the pod strain modes for the hprfe2 model, and the reconstruction of the displacement and damage fields of the rve.
Snapshots in the elastic regime are processed separately from those obtained during the inelastic regime. This strategy guarantees an exact response of the hprfe2 model in the elastic regime.
The sampling was performed using KratosMultiphysics simulation software Ferrándiz et al. (2020).
4.3 Accuracy: assessing nonsampled trajectories
The accuracy of the model is primarily assessed by comparing, for a particular nonsampled validation trajectory, the homogenized stresses obtained with the reduced model and the hf simulation.
The bvp of the microcells of models A and B are solved with the hprfe2 technique using a combination of modes in the range of to , and cubature points in the range of to . The obtained results correspond to an imposed monotonic macrostrain trajectory, defined by
(21) 
until reaching the time . This trajectory is not included in the sampling process.
The error for each of these combinations is defined as
(22) 
i.e., the norm of the maximum relative difference between the homogenized stress tensor given by the reduced model , and the hf model stress , computed along the trajectory. Errors, given in terms of the number of cubature points and strain modes, are plotted in Figures 2 and 3 for the microcells of models A and B. In general, and as expected, we can observe that a higher number of strain modes reduces the approximation error provided that a high enough number of cubature points is used. Also, as expected, the curves tend to an asymptotic value coinciding with the error provided by the standard rom and a full Gauss quadrature rule determined by the fe mesh in the microcells.
Note that by taking a small number of strain modes ( modes) and cubature points (less than ), the assessed errors for both cases (¡]) are much lower than the typical admissible tolerances in engineering problems.
4.4 Model conservation: assessing the physics of the hprfe2 reduced model
The roec scheme used to integrate the balance equation (5) requires the computation of stresses and the damage evolution (6) in the specific points of the microcell, determined by the roec technique. This feature of the hprfe2 model confers a physical basis with an inherent capacity for modeling material responses not specifically considered during the sampling stage, but that still obey basic physical assumptions made at the microscale (e.g., constitutive model thermodynamic consistency, loadingunloading conditions). In this sense, the physics of the reduced model is identical to that shown by the conventional fe approach, where the material physical laws are strictly satisfied in the quadrature points.
To test this feature, the microcells of the models A and B are solved with the hprfe2 technique using strain modes and cubature points. A cyclic macrostrain trajectory is imposed, defined by the macrostrain (21) and the parameter increasing monotonically from to and back to .
Different homogenized stress components are depicted in Figures 3(a) and 4(a) for both microcells. Continuous curves depict solution paths of the hf model while dashed curves correspond to the reduced model. In both cases the solution of the reduced model is remarkably coincident with the hf solution even during the unloading regime, which was not tested in any of the sampling paths commented in the previous section. This result evidences the close connection between the reduced model and the expected physical response of the composite.
Deformed microcells at maximum load are depicted in Figures 3(b) and 4(b). These deformed configurations have been obtained with the reconstructed displacements fields, as indicated in Section 3.3. In the same figures, the reconstructed damage and stress component are plotted on the deformed cells. The reconstructed fields are bounded by an error of 3.24% for model A, and 6.24% for model B, compared with the fields obtained with the hf technique. Note the very high accuracy of the hprfe2 solutions to display the strain jumps caused by the fibermatrix decohesion and interply delamination mechanisms.


4.5 Computational performance
Gains in computing time of hprfe2 with respect to the standard fe2 is quantified by the speedup curves shown in Figure 6.
Speedup curves correspond to 3D microcells B1 to B5, of increasing mesh size with , , , and elements (approximately hexahedra and 6ip wedges). The microcells are based on the described model B, and present a random distribution of fibers (volume fraction of ) with a periodic boundary, generated with the algorithm described in Melro et al. (2008). Each speedup value is computed as the ratio between the time of the hf analysis and the corresponding reduced model analysis, varying in the number of modes and cubature points. hf CPU times are , , , and for models R1 to R5, and CPU times demanded by the reduced model are , , , and for combinations 20m100ip, 30m200ip, 40m400ip, 50m800ip and 60m1600ip, respectively. In the case of reduced analysis, time depends only on the number of modes and cubature points, indistinctly of the RVE that they represent.
As expected, the smaller the number of modes and cubature points, the lower the accuracy and the higher the speedup. Also, the speedup is higher for larger microcells and a given number of strain modes and cubature points. Consequently, in order to determine a fair measure of the speedup for engineering purposes, the maximum admissible error of the evaluated results has to be fixed as a target, which in turn limits the minimum number of strain modes and cubature points that can be taken in the reduced model.
In this view, a series of tests involving microcells of increasing mesh size ( , , , and hexahedra in microcells R1 to R5, respectively) and complexity (defined as the number of algebraic operations) has been previously performed by the authors in LloberasValls et al. (2019), and the resulting speedups are reproduced in Figure 7 (involved CPU times are shown later in section 4.7.2). For every case of the study, the errors obtained by reduced models when taking at least strain modes and cubature points are below (an acceptable error threshold for engineering purposes). Note that periodicity of the solutions is specifically broken by the fact that the geometry is not repetitive, i.e., the larger microcells are not the tiling of the smaller ones. Furthermore, boundary conditions compatible with minimal kinematic constraints are defined, so that the physical complexity increases with the volume (e.g., the strain field of the 4fiber4layer microcell R4 is not the periodic repetition of the strain field of the 2fiber2layer microcell R2.) Thus, the strain modes of the larger microcells are not the periodic repetition the smaller ones.
It is remarkable that, for large microcells with millions of Gauss integration points within the fe discretization, the speedup of the reduced model reaches values in the order of (plus any other speedup gain due to parallel processing). In other words, hf multiscale simulations that would take years of computing time, can now be performed in a few hours with errors below .
4.6 Material customization: dependency on material properties values for the same material model
Several tests are outlined to prove the customization of the methodology in terms of changes in the material parameters characterizing the composite phases. The idea behind these numerical experiments is to evaluate the capability of the hprfe2 technique—specifically developed using a given composite morphology and a given set of material parameters—for simulating composites having the same morphology and ruled by the same constitutive model, but with remarkably different values of material properties with respect to those adopted for the sampling stage. A suitable model response for this kind of test may be a highly desired attribute in applications of material design, as well as in virtual testing.
The microcells employed in the previous sections have been constructed with sampling procedures that have considered the reference material parameters depicted in Tables 1 and 2. In the following numerical experiments, the material response of those reference materials are compared to the same reduced models, but in this case using a different set of material properties without performing new sampling. The microcells of the models A and B are now characterized by the matrix and cohesive interphase parameters displayed in Table 3 and denoted materials Custom M1, Custom M2. Note that Custom M3 is only used in the case of the microcell of the model B, in which the fibers are subjected to damage. In the table, the new material parameters are compared to the reference parameters used in the sampling of the reduced models.
Custom reduced model responses are compared to those obtained using hf models. Result curves versus for each custom material are reported in Figure 8, for hprfe2 bases with strain modes and cubature points. Curves correspond to an imposed cyclic macrostrain trajectory, defined by the macrostrain (21), with increasing monotonically from to (microcell A) and (microcell B), and back to 0. The results obtained with the original material parameters used for the sampling, and for a similar trajectory, are also included in the figure.
These plots show a close agreement between the solutions obtained with both methodologies. A slight difference is observed in Figure 7(a) for microcell of the model A, Custom M2 material. Possibly, this effect is due to the extremely large change in the customized hardening law.
Note that in the case of the curve CM3 in Figure 7(b), the fiber is allowed to damage even though the hprfe2 strain bases used in this analysis has been obtained considering the fibers elastic. Even so, the comparison with the hf is remarkable good.
Original  Custom M1  Custom M2  Custom M3  
Model A  
Matrix:  Young modulus E  5000  –  
Poisson ratio  –  
Elastic threshold  300  –  
Inf elastic threshold  320  –  
Hardening parameters ,  ,  0.5, 0.1  ,  –  
Cohesive,  Young modulus E  –  
Interply:  Poisson ratio  –  
Elastic threshold  –  
Inf elastic threshold  –  
Hardening parameters ,  ,  ,  0.9, 0.9  –  
Model B  
Matrix:  Young modulus E  8000  8000  
Poisson ratio  
Elastic threshold  
Inf elastic threshold  140  
Hardening parameters ,  ,  ,  ,  ,  
Inclusion:  Young modulus E  
Poisson ratio  
Elastic threshold  –  –  –  60  
Inf elastic threshold  –  –  –  70  
Hardening parameters ,  –  –  –  0.01, 0.01 
4.7 Virtual testing of a composite coupon
The coupon depicted in Figure 8(a) is used to perform several virtual tests. It is composed by a unidirectional composite ply, characterized by the microcell B. Three cases are simulated, corresponding to fiber orientations (longitudinal), and (transversal) in the plane , obtained by standard rotation of the microcell to adjust local and global reference axis. The thickness of the specimen is three orders of magnitude larger than the fiber diameter, so the separation of scales required by fe2 method is ensured.
A uniform displacement is imposed on a stiff plate stuck to the inferior surface of the specimen, as shown in Figure 8(b). Exploiting the symmetry of the problem, a fourth part of the coupon is simulated with linear hexahedral elements and Gauss integration points per element. Each integration point is associated with the microcell described previously in Table 2.


The mechanical response of the specimen for different orientations of the fibers is shown in Figure 10. Solutions of reaction forces versus displacement are plotted for three combinations of microstrain modes and cubature points: 20 modes100 cubature point (20m100cp), 30 modes200 cubature points (30m200cp) and 40 modes400 cubature points (40m400cp). Due to the unfeasibility of performing an actual fe2 simulation, the solution taken as reference is the one obtained with the hprfe2 technique, adopting microstrain modes and cubature points.
For the three orientations, the responses of the reduced model agree very well with the reference solutions, with relative errors (with respect to the case 50m1800cp) between and . The larger relative error is observed for the case , 30m200cp. Note that the error, with respect to the reference curve, observed in the case of does not monotonically decrease with the increment of the number of modes and cubature points. According to the conclusions drawn from Figure 3, we argue that the responses of Figure 10 are obtained with a low number of modes and cubature points, in correspondence to a region of nonmonotonic error convergence. However, being that the involved relative errors are small for the purposes of this work, it does not deserve a further analysis involving increments of the number of modes and cubature points.
Figure 11 shows the displacement and stress of the coupon in the mentioned cases at maximum load. It also shows the reconstructed microstress and damage field of the deformed microcell (with the procedure explained in Section 3.3), located at an arbitrary Gauss integration point of the coupon. These microscopic results make evident one of the main features of the fe2 approaches, and preserved in the hprfe2 model, i.e., the detailed analysis the microstructural field, which can not be assessed by simpler phenomenological macroscopic models.



4.7.1 Virtual testing of a composite coupon: customized materials
The previous virtual test is performed with customized material parameters, different from the values used during the sampling stage. Customized parameters correspond to material Custom M2 of Table 3.
The computed results are displayed in the Figure 12 for fiber orientations , and . They are obtained from the hprfe2 model sampled with the original material parameters, but using the custom material parameters during the virtual testing stage. These results are plotted for the cases 20m100cp, 30m100cp and 40m400cp.
The curves denoted as “Reference” are the result obtained with a hprfe2 model sampled with the custom material parameters. They are taken as the reference curves.
To compare the structural effects induced by the change of the customized parameters with respect to the original ones adopted in this test, the curves denoted as “Original” are plotted in the same figure. These curves, which are the reference curves in Figure 10, are obtained with an hprfe2 model, sampled with the original material parameters and using the same material parameters in the virtual testing stage.
For the three orientations, the responses of the reduced model agree very well with the reference solutions, with relative errors between and . The larger relative error is observed for the case , 30m200cp.
4.7.2 Virtual testing of a composite coupon: estimated performance comparison
The same composite coupon shown in Figure 9 is tested using a sequence of microcells employed in the speedup curves in Figure 7. As discussed before in Section 4.5, in terms of mechanical deformation modes, size implies a complexity increment due to their geometry and the use of minimal kinematic boundary conditions.
The sequential computational time for a multiscale analysis of the coupon is reported in Figure 13
, for the estimated time of fe2 and the measured time of the hprfe2 analyses.
For the fe2 cases, time is estimated as the product of the mean iteration time, the mean number of iterations per step, the number of steps and the total number of Gauss integration points of the coupon discretization. The fe2 analysis corresponding to the most complex microcell would demand a year of computation.
The measured time corresponding to the hprfe2 analyses is close to three hours. It should be noted that all reduced bases employed for microcells R1 to R5 involve 40 strain modes and 200 cubature points, yielding relative errors below 1% in the 5 cases. The error analysis is performed in independent analyses considering only one macroscale point.
Consequently, based on this simple multiscale analysis, it is concluded that the hprfe2 presented in this work clearly breaks the barrier imposed by the multiplicative cost within hierarchical multiscale analysis, and represents a real chance to export the fe2 technology to the simulation industry.
5 Conclusions
The hprfe2 technology is assessed in this contribution in terms of its benefits for multiscale solutions for a number of cases of industrial relevance, concerning 3D reinforced composite laminate materials. After an exhaustive evaluation of the presented technique on the above mentioned materials a number of key advantage for the material design industry arise, in comparison to alternative approaches.

Mechanical coherency. The formulation of the multiscale reduced order modelling in terms of the strain field fluctuation turns to be ideal in the context of a fe2like approach, in which downscaling and upscaling contain information of the macrostains and stresses. The strain bases contain only information relevant to the mechanics inheriting the type of boundary conditions employed during the sampling stage, and that is compatible with the strain field, and neglecting the rigid body modes which do not contribute to the constitutive relation. Moreover, the strain snapshots are directly Gauss points quantities and not nodal values, which simplifies the formulation and avoids interpolating quantities to the sampling points needed to integrate the strain modes. Another important feature of the hprfe2 relies in the methodology adopted to integrate the strain modes. The roec integration technique arises from approximating the free energy field (fundamental variational principle) and not its derivatives, as proposed in Hernandez et al. (2017). As a result, the problem goal, i.e., find the minimum of the variational principle (free energy), is the actual key of the problem solution unlike in alternative methods based on finding null values of the functional derivatives. Consequently, it results in a lowcost and very accurate formulation, leading to the high performance reduction technique hprfe2.

Numerical consistency. The accuracy of the approximation, as compared to the hf analysis, improves with the number of the reduced strain modes considered in the construction of the strain basis. Additionally, an increase of integration points (or energy modes), involved in the integration of the strain modes, provides an improvement of the accuracy with respect to the hf solution and, in the limit, tends to the solution provided with a rom model including all quadrature points of the original fe discretization. Noteworthy accuracy is found by adopting a low number of strain modes and cubature points.

Model completeness. Since the material parameters employed within the integration of the reduced strain basis are obtained from the evaluation of the constitutive relation, the link with the physics and the consistent thermodynamics is never lost. For this reason, the reduced model can capture well those loading situations (and even thermodynamic irreversible phenomena, like irreversible loadingunloading processes) that were not included during the sampling stage, without the need for including extra modes, as it would happen in other methodologies (e.g., neural network surrogates).

Material customization. The reduced basis can be seen as a discretization support for the integration of the solution. In this view, one can keep the integration support and substantially change the model parameters for material design purposes. This has been shown in the above sections for significant changes of the material parameters without compromising the quality of the solution. Indeed, this can not be taken as a general recipe and radical changes of the model parameters may, in some cases, demand a more complete strain basis to properly capture the correct solution. However, for the type of materials adopted in this study, the material parameters could be remarkably varied without noticing a clear deviation from the hf solution. The possibility of a reduced model supporting a certain degree of material customizations is considered an overriding quality when the reduced model has to be applied to the design of new materials.

Microstructure monitoring and design. The main feature of fe2 techniques—modeling the coupled physics in two separated scales—is preserved in the presented approach. It relies in monitoring the microstructural behaviour at different key macrostructural points to determine the performance of different phases submitted to stress concentrations and complex loading. This is performed by reconstructing the displacement, stress and internal variable fields as indicated in the formulation of the model. On this view, an “overall” picture of the material behavior at the low scale is obtained, and relevant decisions concerning the microstructural topology and constituents can be taken in order to improve the macrostructural behavior. Again, a desirable feature for the simulation tools used by the material design industry.

High computational performance. A number of tests have been conducted on an engineering coupon with microstructures of increasing increasing size and discretization features. In the presented analyses, the speedup of the approach grows with increasing complexity of the analyzed microstructure, and can reach up to for microstructures with integration points. For the case of fullyperiodic microcells treated with periodic boundary conditions, this result is not really meaningful since volumes are increased following the periodicity in the three spatial directions and, therefore, essentially contain the same mechanical information. However, for the case of minimal kinematic boundary conditions and a sequence of microcells which are not necessarily periodic, a similar observation is done in terms of speedup. It is concluded that for the type of materials adopted in the present study, relevant speedup values can be obtained, which allow computations that would typically span several years to be resolved in a few hours and, therefore, being affordable for industrial purposes.
To the authors’ knowledge, the presented technology represents a novel contribution towards realistic, accurate and affordable industrial multiscale simulations. It promises to advance the current standards for material simulation due to a tangible reduction in experimental testing, more complex material behavior that can be modeled with affordable and reliable numerical simulations, seamlessly integrated in current industrial workflows.
Acknowledgements
The authors acknowledge financial support from the Spanish Ministry of Economy and Competitiveness, through the “Severo Ochoa Programme for Centres of Excellence in R&D” (CEX2018000797S) and the research grant DPI201785521P for the project “Computational design of Acoustic and Mechanical Metamaterials” (METAMAT). This research has also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Proof of Concept Grant agreement 874481) through the project “Computational design and prototyping of acoustic metamaterials for target ambient noise reduction” (METACOUSTIC). The authors also acknowledge the guidance and assistance with the microcells meshes from Dr. Pedro Camanho and Dr. Fermín Otero from INEGI (Portugal) during the preparation of this manuscript.
References
 Variational foundations and generalized unified theory of rvebased multiscale models. Archives of Computational Methods in Engineering 23 (2), pp. 191–253. Cited by: §1.
 High performance reduced order modeling techniques based on optimal energy quadrature: application to geometrically nonlinear multiscale inelastic material modeling. Archives of Computational Methods in Engineering 26 (4), pp. 771–792. Cited by: §1, §1.
 Nonlinear model reduction via discrete empirical interpolation. SIAM Journal on Scientific Computing 32 (5), pp. 2737–2764. Cited by: item iv).
 PGDbased computational vademecum for efficient design, optimization and control. Archives of Computational Methods in Engineering 20 (1), pp. 31–59. External Links: Document, ISBN 18861784 Cited by: item iii).
 A short review on model order reduction based on proper generalized decomposition. Archives of Computational Methods in Engineering 18 (4), pp. 395. Cited by: item iii).
 Symmetric iterative interpolation processes. Constructive Approximation 5 (1), pp. 49–68. External Links: Document, ISBN 14320940 Cited by: item ii).
 A pgdbased multiscale formulation for nonlinear solid mechanics under small deformations. Computer Methods in Applied Mechanics and Engineering 305, pp. 806 – 826. External Links: ISSN 00457825 Cited by: item iii).
 KratosMultiphysics External Links: Document, Link Cited by: §4.2.
 Multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Computer Methods in Applied Mechanics and Engineering 183 (34), pp. 309–330. Cited by: §1.

Accelerating multiscale finite element simulations of historydependent materials using a recurrent neural network
. Computer Methods in Applied Mechanics and Engineering 357, pp. 112594. Cited by: item i).  Representative volume: existence and size determination. Engineering Fracture Mechanics 74, pp. 2518–2534. Cited by: §1.
 Dimensional hyperreduction of nonlinear finite element models via empirical cubature. Computer methods in applied mechanics and engineering 313, pp. 687–722. Cited by: item iv), §1, §3.2, item 1.
 Highperformance model reduction techniques in computational multiscale homogenization. Computer Methods in Applied Mechanics and Engineering 276, pp. 149–189. Cited by: item iv).
 A multiscale method for periodic structures using domain decomposition and ecmhyperreduction. Computer Methods in Applied Mechanics and Engineering 368, pp. 113192. Cited by: item iv).
 Generation of energyminimizing point sets on spheres and their application in meshfree interpolation and differentiation. Advances in Computational Mathematics 45, pp. 3021–3056. Cited by: §4.2.
 The latin multiscale computational method and the proper generalized decomposition. Computer Methods in Applied Mechanics and Engineering 199 (2122), pp. 1287–1296. Cited by: item iii).
 Reduced finite element square techniques (rfe2): towards industrial multiscale fe software. In Proceedings of the XV International Conference on Computational Plasticity  Fundamentals and Applications, Cited by: §1, §1, §4.5.
 A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. Journal of Computational Physics 330, pp. 192–220. Cited by: §1.
 Generation of random distribution of fibres in longfibre reinforced composites. Composites Science and Technology 68 (9), pp. 2092–2102. Cited by: §4.5.
 Effective properties of composite materials with periodic microstructure: a computational approach. Computer methods in applied mechanics and engineering 172 (14), pp. 109–143. Cited by: §1.
 Computational micromacro transitions and overall moduli in the analysis of polycrystals at large strains. Computational Materials Science 6, pp. 372–382. Cited by: §1.
 On the existence of representative volumes for softening quasibrittle materials  a failure zone averaging scheme. Comput. Meth. App. Mech. Eng. 199, pp. 3028–3038. Cited by: §1.
 Reduced order modeling strategies for computational multiscale fracture. Computer Methods in Applied Mechanics and Engineering 313, pp. 560 – 595. Cited by: §1, §1, §2, §3.2, §4.2.
 Continuum approach to computational multiscale modeling of propagating fracture. Computer Methods in Applied Mechanics and Engineering 294, pp. 384–427. Cited by: Appendix A.
 Micromechanicsbased surrogate models for the response of composites: a critical comparison between a classical mesoscale constitutive model, hyperreduction and neural networks. European Journal of Mechanics  A/Solids 82, pp. 103995. External Links: ISSN 09977538 Cited by: item i).
 Hyperreduction of mechanical models involving internal variables. International Journal for Numerical Methods in Engineering 77 (1), pp. 75–89. Cited by: item iv).
 A class of general algorithms for multiscale analyses of heterogeneous media. Computer methods in applied mechanics and engineering 190 (40), pp. 5427–5464. Cited by: §1.
 Wavelet based reduced order models for microstructural analyses. Computational Mechanics 63 (3), pp. 535–554 (English). External Links: Document, ISSN 01787675 Cited by: item ii).
 A multilevel projectionbased model order reduction framework for nonlinear dynamic multiscale problems in structural and solid mechanics. International Journal for Numerical Methods in Engineering 112 (8), pp. 855–881. Cited by: item iv).
Appendix A Damage model
A damage model with an elastic domain accounting for tractiononly and a bilinear hardening law is adopted for describing the constitutive response of the components at the microscale. The main features of the model are described in Table 4. The employed constitutive model is extensively described in Oliver et al. (2015) for the case of linear softening behaviour. However, a short explanation of the parameters and symbols used in the table is included in the following lines: is the free energy and is the Hooke’s elasticity tensor written in terms of the Lamé parameters and , with and the forth and secondorder identity tensor. The isotropic damage variable depends on the internal variable . This internal variable determines the size of the elastic domain through the damage function . The term in the definition of is a generalized norm of the strains accounting for the positive effective principal stresses . A bilinear hardening law is considered in the present study with two values of the hardening parameter . It is worth mentioning that no regularization is required since softening behaviour is not included in the present study.
Free energy :  

Damage variable :  
Constitutive equation:  
Damage function :  
Initial conditions:  
Loadingunloading conditions:  
Hardening law : 