* At the macroscale level, were adapted the last advances of the ''Continuum Strong Discontinuity Approach (CSDA)'', developed for monoscale models, devising a new finite element exhibiting  good ability to capture and model strain localization in bands which can be intersect the finite element in random directions; for failure propagation purposes, the adapted ''Crack-path  field'' technique oliver2014crack, was used.

First, the standard microscale finite element model ''High Fidelity (HF)'', is projected and solved in a low-dimensional space via ''Proper Orthogonal Decomposition (POD)''. Second, two  techniques have been developed and studied for multiscale models, namely: a) ''interpolation'' methods, and b) ''Reduced Order Cubature (ROQ)'' methods An_2009. The reduced bases for the projection of the ''primal variables'', are computed by means of a judiciously training, defining a set of pre-defined ''training trajectories''.

BVPBoundary Value Problem        CSDAContinuum Strong Discontinuity Approach        DNSDirect Numerical Simulation        EBAExpanded Basis Approach        EFEMEmbedded Finite Element Methodology        EIMEmpirical Interpolation Method        FEFinite Element Method                  FE<math>^2</math>FE<math>\times </math>FE hierarchical multiscale technique        HFHigh Fidelity        HPROMHigh-Performance Reduced Order Model        HROMHyper-Reduced Order Model        MORModel Order Reduction        PODProper Orthogonal Decomposition        ROMReduced Order Model        ROQReduced Optimal Quadrature        RUCRepeating Unit Cell        RVERepresentative Volume Element        SVDSingular Value Decomposition

In this section, a general insight about multiscale approaches, failure modeling and Model Order Reduction, is given. The aim here is to describe the most important contributions along the  history of the continuum mechanics community in these research fields.

During the last decades, a large variety of multiscale strategies focusing on the study and analysis  of the mechanical behavior of heterogeneous materials, have been proposed by the computational mechanics community. Based on the work of Bohm_2013,  these strategies may be divided into three main groups:

* ''Mean-Field approaches'': based on the seminal contributions developed by Eshelby_1957 and Mori_Tanaka_1973. In these approaches, the microfields within each constituent of  an heterogeneous material, are approximated by their phase averages, typically, ''phase-wise'' uniform stress and strain fields are employed. Recently, the application of these approaches to nonlinear  modeling in composites has become a subject of active research.

The formal treatments were provided by, ,  Nemat_Nasser,Bornert_1996,Ponte_Suquet_1998,Markov_2000,Milton_2002,Torquato_2002. Two of the most relevant results of this kind of models are,  the ''upper bounds'' of Voigt_1889, and the ''lower bounds'' of Reuss_1929. Posteriorly, Hill extends those bounds to tensorial entities, particularly, for constitutive tangent tensors Hill_1952.

Generally speaking, these bounds, while universal and very simple to assess, offer poor approximations. The reason relies on the fact that these measures do not contain any information on the geometry of the different phases beyond the phase volume fraction, being too slack for many practical purposes. However, considerably tighter bounds on the macroscopic behavior can be obtained from a variational formulation Hashin_Shtrikman_1963.

In consequence, it can be built a parametric grid based on the strain space whose nodes are  related with homogenized stress tensors produced by the microscale mechanical response; in order to obtain a random solution, these nodal values can be interpolated via Neural Networks, or Finite Element approximations whose interpolation level are based on the spatial variation of the homogenized stress tensor. In recent works, , Ferrer_et_Al_2016, this methodology is applied to multiscale shape optimization, in which, the set of precomputed strain histories, forming a parametrized spherical domain, are associated to their corresponding homogenized stress states, and  microscale optimal shapes.

However, the effectivity of this methodologies is based on the dimensionality of the input parameters, the sampling becomes complex for high-dimensional input parameters, then, those methodologies are inappropriate for microscales involving complex morphologies and/or complex phenomena. Apart from the previous application in multiscale optimization, the most representative methodologies are  presented in Gurson_1975,Tvergaard_1981,Giusti_Blanco_Gurson

* ''Hierarchical models:'' the RVE, subjected to consistent boundary conditions<span id="fnc-1"></span>[[#fn-1|¹]] is used to obtain a detailed microscale response. The link between scales is reached by means of an energetic identity, such as the Hill-Mandel Principle of Macro-Homogeneity Hill_1965,Mandel:1971, or even on more general approaches, like the Principle of Multiscale Virtual Power Blanco2016. In the context of the Finite Element Method  (FEM), this methodology is known as FE<math>^2</math> Feyel_Chaboche_2000.

In virtue of the potential applications in microstructures with complex morphologies, the FE<math>^2</math> technique is suitable to deal with problems facing material instabilities, like fracture processes. In  this sense, some approaches have been proposed Belytschko_multiscale1,Belytschko_multiscale2,Nguyen_et_al_2010,Nguyen_et_al_2011,Matous_et_al_2008,coenen2012multi,Toro_al_FOMF_2014, among  others. One of the main contributions of this work, is an alternative approach with marked differences with respect to the previous ones.

Like the above mentioned models, concurrent models are widely applied. In approaches based on the Finite Element Method, the embedded meshes are not necessarily compatible. However, its computational  viability is only for cases with small scale separation,  this feature becomes into its main disadvantage. Some concurrent approaches in the field of fracture mechanics have been recently developed, , lloberas2012multiscale.

Procedures based on the concept of free energy anderson, becomes into one of the early works in the field of Fracture Mechanics. The concept of ''fracture energy'' appears as a  consequence of this approach.

* ''Cohesive models:'' Based on the introduction of interfaces embedded into a continuum medium. These interfaces admit the development of displacements discontinuities. Cohesive forces  across the interfaces act opposing to the crack opening, diminishing as material degradation takes place. The energy necessary to produce a crack is equal to the fracture energy. Some applications  of this kind of models are found in pandolfi1999finite,molinari2007cohesive,Toro_al_2014,rodrigues20162d.

* ''Continuum regularized models:'' Characterized by a continuum constitutive law displaying a softening response. These kind of approaches are subjected to material instabilities and bifurcation processes, causing ill-possedness of the problem from the mathematical point of view. As a remedy to this flaw, a constitutive law regularization (localization  limiters) is introduced, ensuring mesh objective solutions. Some proposed models can be found in  Pijaudier_Cabot_Bazant_1987,Tvergaard_et_al_1995,Pijaudier_Cabot_Bazant_1987,Aifantis_1984b,Borst_Muhlhaus_1992,Peerlings_et_al_1996,Peerlings_et_al_2001,Peerlings_et_al_2002,Steinmmann_et_al_1991,Muhlhaus_et_al_1987, among others.

* ''CSDA:'' This approach establishes a link between ''cohesive models'' and ''continuum models''. Its fundamentals have been presented in the seminal work Simo_et_al_1993,  posteriorly improved and applied to many applications in static and dynamic scenarios Oliver_2000,Oliver_et_al_2002,Oliver_Huespe_2004a,Oliver_Huespe_2004b,oliver2014crack.  This approach provides an unified theory, which goes from the continuum description to the degradation and posterior material failure exhibiting displacement discontinuities (cracks). In this  context, the continuum constitutive model subjected to a kinematics inducing displacement discontinuities represents also a “projected” cohesive law on the crack surface.

The study of heterogeneous materials subjected to softening, and, therefore, to degradation and failure, through multiscale approaches brings additional challenges. The fundamental reason lies in  two aspects: (a) it becomes imperative the use of regularized constitutive theories at both scales in order to ensure the well-possednes of the multiscale problem. (b) The ''size effect'',  intrinsically related to the fracture energy, and extensively studied by Bazant_Planas_1998. As a result of this, the homogenized stress tensor, in the post-citrical regime,  becomes extremely sensible to the RVE size.

The second issue is the existence of the RVE Gitman_et_al_2007, and the fulfillment of the basic hypothesis in multiscale modeling.

The necessity to develop specific homogenization techniques, becomes a starting point for obtaining consistent multiscale formulations. Belytschko_multiscale1 has proposed a methodology that  excludes the localization domain in the homogenization process. More recently, Belytschko_multiscale2 proposed a predetermined size of RVEs. Matous_et_al_2008 describes a novel methodology, based on the existence of a macroscopic adhesive interface, which links the macroscale jump of displacements with an equivalent jump at the microscale, imposed by  consistent boundary conditions.

Recent works Nguyen_et_al_2010,Nguyen_et_al_2011, describe the material failure by means of nonlocal gradient theories. In this kind of approaches,  and, in contrast with other alternatives,  the homogenization of the stress field during the post-critical regime, is carried out at the localization zone (this zone corresponds to a subdomain of the RVE). However, other authors claimed  some inconsitencies related to this kind of approaches, particularly, about the fact that kinematics at the macroscale is not equivalent to the kinematics modeled at the microscale.

Toro_al_FOMF_2014 presents a more consistent theory which some several shortcomings given by the previous approaches, however, this the homogenization procedure evolves along time, producing an algorithm too invasive.

Although no doubt the most versatile and accurate homogenization technique, with no other limitation in scope than the hypotheses of statistical homogeneity and  scale separation, the direct computational homogenization approach violates squarely the modeling precept outlined at the outset &#8211;it does not discriminate between essential and irrelevant features in solving the fine-scale ''BVP''&#8211;, making the accuracy/parsimony balance to tilt unduly towards the accuracy side and far from the parsimony one. The consequence is its enormous computational cost (in comparison with analytical and semi-analytical homogenization techniques).

The idea of exploiting the combination of dimensionality reduction and multiscale modeling is certainly not new. A survey of the related literature reveals that, over the last decade, researchers from various scientific disciplines dealing with multiscale problems have begun to consider the model reduction as a potential route &#8211;complementary to improvements in software and hardware power &#8211;to diminish the often unaffordable cost of multiscale simulations. In the specific context of homogenization-based multiscale methods, the application of model  reduction techniques has been addressed by several authors, namely, Ganapathysubramanian_2004,yvonnet2007reduced,Boyaval_2007,Monteiro2008,Nguyen_ROM_2008. The strategy adopted in all  these works for constructing a cost efficient model of the micro-cell is the standard reduced basis method. The gist of this strategy is to project the governing equations onto a low-order subspace  spanned by carefully chosen bases Amsallem_2009.

The origin of the first effective proposal on this issue can be traced back to the seminal work of Barrault_2004, who suggested to approximate the nonlinear term in the reduced-order equations by a linear combination of a few, carefully chosen basis functions. In the spirit of a offline/online strategy, in the standard reduced basis approach, these spatial bases are computed offline from full-order snapshots of the non-linear term, whereas the corresponding parameter-dependent modal coefficients are determined online by interpolation at a few (as many as basis functions),  judiciously pre-selected spatial points. As in classical reduced bases methods, the efficiency of this second or collateral reduction is predicated on the existence of a moderate number <math display="inline">M \ll N</math> (<math display="inline">N</math> is the original dimension of the problem) of basis functions whose span accurately approximate the manifold induced by the parametric dependence of the nonlinear contribution.  The interpolation method developed by Barrault_2004 is known as the EIM; the main ingredients of this method are: a) the use of a greedy algorithm to  generate a set of maximally independent bases from the collection of snapshots of the nonlinear term, on the one hand; and b) the recursive selection &#8211; also via a greedy algorithm &#8211; of  spatial locations where the error between the full-order bases and their reconstructed counterparts is maximum<span id="fnc-2"></span>[[#fn-2|¹]].

In solution methods in which the governing equations are used in its variational form (as in the FE), the reduction of the computational complexity arising from nonlinearities can  be, alternatively, achieved by approximating the integrals in which the offending nonlinear function appears, rather than the function itself, as done in the interpolatory and least-square  reconstruction techniques discussed above. Based on this observation, An_2009 propose a quadrature scheme devised for fast-run integration of the subspace spanned by a representative  set of snapshots of the nonlinear integrand.

The development of reduced models for non-homogeneous materials has been tackled in numerous previous contributions, such as michel2001computational, where the proposed reduction techniques are based on Fourier's transforms, or yvonnet2007reduced, where a reduced model is applied the homogenization analysis of hyperelastic solids subjected to finite strains. Also, the work in Ryckelynck2009hyper develops a hyper-reduced model of a monoscale analysis which consider nonlinear material behavior. However, the existing literature barely considers reduced order modeling of non-smooth problems, as is the case of fracture, where discontinuous displacements occur. The multiscale case, when fracture takes place at both scales of the problem, makes the task even  much harder. Indeed, only very few contributions have been presented in the literature about this topic, see for example: oskay2007eigendeformation, which follows an eigendeformation-based  methodology, or zhang2016reduced,kerfriden2013partitioned that resort to global&#8211;local approaches.

As pointed out earlier, many ways of approaching the study of microscale behavior has been proposed. Particularly, two main procedures dominate the scenario of hierarchical multiscale modeling at  present: analytical approaches mainly based on the ''asymptotic expansion analysis'',  and numerical approaches based on variational formulations supported by physical arguments. The later relies on the concept of ''Representative Volume Element (RVE)'', being the ''Computational Homogenization'' approach, the most common one. The model reduction techniques developed in this work will be applied to hierarchical multiscale models derived from the latter approach. Their fundamental assumptions are presented in its more conventional form.  For a more in-depth description of the underlying axiomatic framework, the readers are referred to Miehe_et_al_1999,Feyel_Chaboche_2000,Kouznetsova_thesis,Feijoo_multi_scale_2006,DragoPindera:2007,WeinanE:2011

This representative sample, hereafter denoted <math display="inline">{\mathcal{B}_{\mu }}\in {\mathbb{R}^d}(d=2,3)</math>, is assumed to exhibit several features. One of those corresponds to the  ''size indifference'' property Terada2000,Kouznetsova_thesis,DragoPindera:2007, which states that if the size of this sample is increased, the response remains identical regardless  the ''admissible'' boundary conditions on the RVE. The lower size limit for the RVE satisfying the size indifference property is represented by the  ''characteristic length-scale'' denoted as <math display="inline">{h_\mu }</math>, giving rise to the existence of the RVE, whereas in microstructures that display  periodicity, is known as ''RUC'', or simply ''unit cell''. Furthermore, <math display="inline">{\mathcal{B}_{\mu }}</math> has to be small enough to be regarded as a point at the macroscale GrossSeelig:2011 (, <math display="inline">{h_\mu }\ll {L}</math>, being <math display="inline">{L}</math> the characteristic length of the macroscale <math display="inline">{\mathcal{B}}</math>, see Fig. [[#img-1|1]]) this is the so-called ''scale separation'' hypothesis.

* The ''Hill-Mandel Variational Principle of Macro-Homogeneity'' Hill_1965,Mandel:1971, which states an equivalence between the virtual power densities between micro- and macroscale,  and requires the adoption of specific kinematically admissible displacement fluctuations at  the RVE, is adopted.

|[[Image:draft_Samper_355780031-MultiscaleSketch.png|390px|Macrosctructure with an embedded local microstructure.]]

<span style="text-align: center; font-size: 75%;">([[#fnc-3|¹]]) However, in one article supporting this work, dynamic problems are also considered. See Hernandez_HPROM_2017</span>

Eq. [[#eq-2.4|2.4]], can be cast four well-known classes of multiscale constitutive  models emphasizing that the different classes of multi-scale models differ from one another only in the definition of the subspace <math display="inline">\mathcal{\tilde V}_{\mu }\subset \mathcal{\tilde V}_{\mu }^{u}</math> where <math display="inline">{\boldsymbol{\dot{{\tilde u}}}_\mu }</math> belongs  Feijoo_multi_scale_2006:

The scale bridging equations are completed by introducing the ''Hill-Mandel Principle of Macro-Homogeneity'' Hill_1965,Mandel:1971. Based on physical arguments, this Principle states that the  macroscopic stress power equates the volume average over the RVE of the microscopic stress power, making both, macroscale and microscale, continuum descriptors ''energetically  equivalent''. Thus, departing from:

Consider a time discretization of the interval of interest <math display="inline">[t_{0},t_{f}]=\bigcup _{i=1}^{\hbox{nstp}}[t_{n},t_{n+1}]</math>. The current value of the microscopic stress tensor <math display="inline">{\boldsymbol{\sigma }}_{\mu }</math> at each <math display="inline">\boldsymbol{y}\in {\mathcal{B}_{\mu }}</math> is presumed to be entirely determined by, on the one hand, the current value of the microscopic strain tensor <math display="inline">{\boldsymbol{\varepsilon }_\mu }_{n+1}={\boldsymbol{\varepsilon }}_{n+1}+{{\boldsymbol{\nabla }^s_{\boldsymbol{y}}}{\boldsymbol{\tilde{u}}_\mu }}</math>, and, on the other hand, a set of microscopic internal variables <math display="inline">{\boldsymbol{\mu }}_{n+1}</math> that encapsulate the history of microscopic deformations.The (incremental) RVE equilibrium problem at time <math display="inline">t_{n+1}</math> can be stated as follows: given the ''initial data'' <math display="inline">\{ {\boldsymbol{\tilde{u}}_\mu }_{n}(\boldsymbol{y}),{\boldsymbol{\varepsilon }}_{n},{\boldsymbol{\mu }}_{n}(\boldsymbol{y})\} </math> and the ''prescribed'' macroscopic strain tensor <math display="inline">{\boldsymbol{\varepsilon }}_{n+1}</math>,  find <math display="inline">{\boldsymbol{\tilde{u}}_\mu }_{n+1} \in \mathcal{\tilde V}_{\mu }^{u}</math> such that HdezEtAlMon:2012:

also the constitutive homogenized tangent tensor <math display="inline">\mathbb{C}_{n+1}^h</math> which is composed by two parts, on the one hand, the contribution of the constitutive tensor <math display="inline">\mathbb{C}^h</math>, and on the other hand,  the contribution given by the fluctuation displacement field <math display="inline">\tilde{\mathbb{C}}^h</math> Feijoo_multi_scale_2006, for the multiscale modeling of propagating fracture is necessary to upscale additional  tensorial variables in order to guarantee the proper dissipation at both scales, and also fulfill objectivity conditions widely known in fracture mechanics, these variables will be detailed in the Section [[#2.2 Multiscale Fracture Mechanics issues|2.2]].

Let <math display="inline">{\mathcal{B}_{\mu }}=\bigcup _{i=1}^{\hbox{n}_e}{\mathcal{B}}_{\mu }^e</math> be a finite element discretization of the RVE, and let <math display="inline">\{ N_1(\boldsymbol{y}) \cdots N_n(\boldsymbol{y})\} </math> (<math display="inline">n</math> denotes the number of nodes  of the discretization) be a set of ''shape functions'' associated to this discretization. Now we approximate <math display="inline">{\boldsymbol{\tilde{u}}_\mu }\in \mathcal{\tilde V}_{\mu }^{u}</math> and <math display="inline">{\boldsymbol{\eta }}\in \mathcal{\tilde V}_{\mu }^{u}</math> as Hughes_1987:

Now, introducing Voigt's notation, and numerical integration via Gaussian quadrature zienk_taylor_2000,Belytschko_WinKL_2001, expression [[#eq-2.17|2.17]] can be expressed  in matrix format as:

* The existence of the RVE has been questioned in the literature, arguing that for fracture cases, the material loses its statistical homogeneity Gitman_et_al_2007,Nguyen_et_al_2010.

* The fact that the homogenized constitutive model lacks an internal length Bazant_2010, raising similar issues than in classical phenomenological monoscale problems.

Additionally, mesh-bias dependence, and the proper fracture energy dissipation issues Rots_1988 via regularized constitutive models  Oliver_1989,Oliver_et_al_2002,Oliver_Huespe_2004b,oliver2015continuum are also crucial issues to be considered at each scale.

At this point, two main aspects have to be dealed in multiscale hierarchical approaches, the possible non-existence of the RVE combined with the lack of objectivity of its representative size.  Nowadays, several studies have been developed in order to analyze the influence of the use of an unit cell rather than the RVE, by assessing the amount of deviation of the apparent properties  obtained by the unit cell modeling Ghosh_1996,Moulinec_Suquet_1998,Michel_et_al_1999. In this work, all those previous issues have been studied and tested. The use of a complex  ''failure cell'' (representing all features of the micro/meso scale) plays an important role for allowing crack nucleation and posterior coalescence. Along this section, the most important  issues of this multiscale approach are summarized, being the analytical details referenced to the corresponding article.

Along this section, the most important aspects of the proposed multiscale approach are summarized. This multiscale approach is fully detailed in Paper <math display="inline">\#{2}</math> in Sec. [[#5.2 Paper <math>\#{2}</math>|5.2]].

|[[Image:draft_Samper_355780031-figure1.png|480px|Macroscopic (Structural scale) body \mathcalB (a) subdivision in a non-smooth domain \mathcalB_\hboxloc(t),  and a smooth domain \mathcalB⧹\mathcalB_\hboxloc(t) (b) ''h''-regularized displacement and strain discontinuity kinematics.]]

{\dot{\boldsymbol{\varepsilon }}}(\boldsymbol{x},t) = \left({\boldsymbol{\nabla }_{\boldsymbol{x}}}\otimes {\dot{\boldsymbol{u}}}(\boldsymbol{x},t) \right)^{s} \equiv {\boldsymbol{\nabla }^s_{\boldsymbol{x}}}{\dot{\boldsymbol{u}}}(\boldsymbol{x},t)  \quad \quad \forall \boldsymbol{x}\in {\mathcal{B}}\setminus {\mathcal{B}_{\hbox{loc}}}(t)  

 

</math>

* Domain <math display="inline">{\mathcal{B}_{\hbox{loc}}}(t)</math>: the set of points exhibiting material failure and, therefore, a non-smooth behavior. See: <p>'''Sec.2.1 </p><p>Paper <math>\# </math>2''' The strain field at these points is  assumed to be captured by a ''h''-regularized ''strong-weak'' discontinuity kinematics, ''h'' being the width of the corresponding strain localization band (see Fig.  [[#img-2|2]]-(b)).

{\dot{\boldsymbol{\varepsilon }}}(\boldsymbol{x},t) = {\dot{\bar{\boldsymbol{\varepsilon }}}}(\boldsymbol{x},t) + \frac{{\kappa _{{\mathcal{B}_{\hbox{loc}}}}}(\boldsymbol{x})}{h}{\dot{\boldsymbol{\gamma }}}(\boldsymbol{x},t) = {\dot{\bar{\boldsymbol{\varepsilon }}}}(\boldsymbol{x},t)+{\delta ^{h}_{s}}(\boldsymbol{x}){\dot{\boldsymbol{\gamma }}}(\boldsymbol{x},t)  \quad \forall \boldsymbol{x}\in {\mathcal{B}_{\hbox{loc}}}(t)  

 

</math>

|[[Image:draft_Samper_355780031-figure1A.png|240px|Colocation function κ_{\mathcalB_\hboxloc}(x)]]

|[[Image:draft_Samper_355780031-figure2.png|600px|Outline of the multiscale model for propagating fracture: '''a)''' macro and micro scales;  '''b)''' microcell model accounting for material failure.]]

{\dot{{\boldsymbol{\sigma }}}_{\mu }}= \Sigma ^{\hbox{hard}}({\dot{\boldsymbol{\varepsilon }}_\mu }) \equiv \mathbb{C}^{\hbox{hard}}_{\mu }:{\dot{\boldsymbol{\varepsilon }}_\mu }\quad \quad \forall \boldsymbol{x}\in {\mathcal{B}}; \quad \forall \boldsymbol{y}\in {\mathcal{B}_{\mu }}\setminus {\mathcal{B}_{\mu ,\hbox{coh}}}   

 

</math>

* The failure cell, <math display="inline">{\mathcal{B}_{\mu }}</math> is associated to a non-smooth material point at the macroscale (<math display="inline">\boldsymbol{x}\in {\mathcal{B}_{\hbox{loc}}}</math>). Without limiting the use of any other inelastic constitutive model, this domain is    endowed with an ''isotropic continuum damage model'', exhibiting inelasticity with regularized strain softening only for tensile stress - ''tensile-damage continuum damage model''    Oliver_1995b,Faria_1998,Oliver_et_al_implex_2006. Its constitutive response isSee: <p>'''Sec.2.2 </p><p>Paper <math>\# </math>2''' represented in a general form as:

{\dot{{\boldsymbol{\sigma }}}_{\mu }}= \Sigma ^{\hbox{inelas}}({\dot{\boldsymbol{\varepsilon }}_\mu },\boldsymbol{\mu }) \quad \quad \forall \boldsymbol{x}\in {\mathcal{B}_{\hbox{loc}}}(t); \quad \forall \boldsymbol{y}\in {\mathcal{B}_{\mu ,\hbox{coh}}}   

 

</math>

{\dot{{\boldsymbol{\sigma }}}_{\mu }}= \Sigma ^{\hbox{elastic}}_{\hbox{inst}}({\dot{\boldsymbol{\varepsilon }}_\mu }) = \mathbb{C}^{\hbox{elastic}}_{\hbox{inst}}:{\dot{\boldsymbol{\varepsilon }}_\mu }\quad \forall \boldsymbol{x}\in {\mathcal{B}}\setminus {\mathcal{B}_{\hbox{loc}}}(t); \quad \forall \boldsymbol{y}\in {\mathcal{B}_{\mu ,\hbox{coh}}}   

 

</math>

Now, in virtue of [[#eq-2.1|2.1]], and for the sake of generality, the microscopic displacement can be modified for multiscale complex problems, which involves propagation of fracture or other complex phenomena at the large, in the following way:

|[[Image:draft_Samper_355780031-figure3.png|270px|]]

|[[Image:draft_Samper_355780031-figure3A.png|270px|Cohesive Band behavior.]]

Eq. [[#eq-2.26|2.26]] constitutes the displacement counterpart of a ''k''-regularized strong discontinuity kinematics Oliver_1996a, and proves that the herein proposed cohesive-bands approach,  is consistent with a ''k''-regularized strong discontinuity at the cohesive domain <math display="inline">{\mathcal{B}_{\mu ,\hbox{coh}}}</math>. In  consequence, the corresponding microscopic strain fluctuation field is given by:

From Eq. [[#eq-2.28|2.28]], it can be concluded, that the second term at the right-hand side becomes unbounded in the limit <math display="inline">k \to 0</math>. In multiscale modeling, this expression is equivalent  to the one given, in phenomenological monoscale models, in the context of the ''Continuum Strong Discontinuity Approach'' (CSDA) of material failure Oliver_et_al_2002.

One of the most specific features of the proposed multiscale approach, is that the same homogenization setting is used in points of both domains, smooth (<math display="inline">\boldsymbol{x}\in {\mathcal{B}_{\mu }}\setminus {\mathcal{B}_{\mu ,\hbox{coh}}}</math>),  and non-smooth (<math display="inline">\boldsymbol{x}\in {\mathcal{B}_{\mu ,\hbox{coh}}}</math>), coinciding with the approach presented in Sec. [[#2.1 Computational Homogenization|2.1]]. Other approaches Toro_al_FOMF_2014, redefine the failure cell along time, fulfilling conditions of material bifurcation induced by instabilities at the microscale. More complex approaches Kouznetsova_thesis,Geers_et_al_2010,Otero_Martinez_2015,Lesicar_et_al_2015 propose the  use of ''second-order'' computational homogenization schemes in order to get better accuracy in the prediction of high strain gradients. In this work it is claimed the ability of the proposed  approach to induce discrete failure in a ''first-order homogenization'' setting, giving rise to objective responses, and proper energy transfer through scales.

An issue appearing in this scenario, widely known in hierarchical multiscale approaches, is its high computational cost. In this context, the proposed model was also conceived to be combined with the use of ''model  order reduction techniques'' (Paper <math display="inline">\#{3}</math>) Oliver_Caicedo_HROM_2017. These techniques have been deeply studied in this work, and their main features are presented in Chapter [[#3 Model Order Reduction in Multiscale Analysis|3]].

|[[Image:draft_Samper_355780031-figure4.png|540px|Multiscale model: '''(a)''' failure cell with activated failure mode; '''(b)''' geometrical characterization of the failure mode.]]

For the sake of generality, the RVE is considered composed by several components: a matrix, and randomly distributed inclusions and voids. In addition, a  number of cohesive bands are considered defining the set <math display="inline">{\mathcal{B}_{\mu ,\hbox{coh}}}</math> (a sketch is presented in Fig. [[#img-6|6]]); those cohesive bands allow failure within the matrix,  See:

 

Paper <math>\# </math>2''' across the aggregates and at the matrix/aggregate interface.

 

Among them, there are <math display="inline">n_{active}</math> cohesive bands, defined as <math display="inline">{\mathcal{B}^{(i)}_{\mu ,\hbox{act}}}</math>, <math display="inline">i=1, \cdots n_{active}</math> which are in an inelastic softening state,  defining a specific ''failure mode'' <math display="inline">{\mathcal{B}_{\mu ,\hbox{act}}}\subset {\mathcal{B}_{\mu ,\hbox{coh}}}</math>, the remaining set of bands are in unloading state, and are discarded for the purposes of this analysis. As a consequence of the  previous statement, the discrete version for the middle-line of the cohesive band <math display="inline">(i)</math> (denoted as <math display="inline">S_\mu ^{(i)}</math>), and the domain defining the failure mode <math display="inline">{\mathcal{B}_{\mu ,\hbox{act}}}</math> can be expressed as follows:

where <math display="inline">S^{(i)}_\mu </math> can be regarded as the middle-line of the cohesive band segment <math display="inline">(i)</math>.

being, <math display="inline">L_{\mu }=\left|S_\mu \right|</math> the measure (length in 2D and area in 3D) of <math display="inline">S_\mu </math> and <math display="inline">\overline{(\cdot )}_{S_\mu }</math> stand for the average value of <math display="inline">\overline{(\cdot )}</math> in the activated microscopic failure mechanism. oliver2015continuum. The term <math display="inline">{\dot{\boldsymbol{\chi }}}(\boldsymbol{x},t)</math> vanishes, in case of components with equal material properties at the microscale, this is the trivial one also known as  ''homogeneous'' case, in addition, depending on the boundary conditions, the displacement fluctuation field is also close to zero.

Finally, after some manipulations of Eq. [[#eq-2.30|2.30]], and following the definitions of microscale kinematics in Eq. [[#eq-2.28|2.28]], and the lemma in Eq. 23 in  Paper <math display="inline">\#{2}</math> oliver2015continuum, the resulting homogenized constitutive equation fulfills the following:

where, <math display="inline">{l_\mu }</math> stands for a ''characteristic length'', depending on the activated microscopic failure pattern. The tensorial entities <math display="inline">{\dot{\boldsymbol{\chi }}}(\boldsymbol{x},t)</math> and <math display="inline">{\dot{\boldsymbol{\varepsilon }}}^{(i)}(\boldsymbol{x},t)</math>, are inelastic strains, and play the same role than internal variables in phenomenological models. However, unlike them, here, See:

 

'''Sec.2.4  

Paper <math>\# </math>2'''  their evolution is determined, at every macroscopic sampling point <math display="inline">\boldsymbol{x}</math>, by homogenized values of entities at the corresponding microscopic failure cell <math display="inline">{\mathcal{B}_{\mu }}</math>. This extends to non-smooth problems, some theoretical results already derived for smooth problems, see Michel_2003,Michel_2004. In addition, a characteristic length <math display="inline">{l_\mu }</math> emerges naturally in Eq. [[#eq-2.36|2.36]],  as the ratio between the measure of the failure cell (area in 2D and volume in 3D), and the measure (length/surface) of the activated microscopic failure mechanism. In  consequence this length is of the order of the failure cell size. For a deeper review of the analytical results of this induced homogenized constitutive model, the reader is addressed to Sec. 2.4 in Paper <math display="inline">\#{2}</math>.

The role of the characteristic length, <math display="inline">{l_\mu }</math>, ''naturally derived'' from the present formulation, is not only computational, but it has also other very relevant physical and mechanical implications. Consideration of such a characteristic length, for multiscale based approaches, has been claimed from the material mechanics community Bazant_2010, and sometimes introduced  in a heuristic way in other approaches Unger_2013. This characteristic length depends on both the specific data of the problem and the local microscopic failure state. Through  its consideration, the correct energy transfer between scales and mesh size objectivity can be achieved.

\overline{(\dot{\boldsymbol{\beta }}_{\mu }\otimes \boldsymbol{n}_{\mu })^s}_{S_\mu } = (\dot{\boldsymbol{\beta }}_{\mu }(\boldsymbol{x},t) \otimes  \overline{\boldsymbol{n}_{\mu }(\boldsymbol{x},\boldsymbol{y})}_{S_\mu })^s = (\dot{\boldsymbol{\beta }}_{\mu }(\boldsymbol{x},t) \otimes \boldsymbol{a}(\boldsymbol{x}))^s     

 

</math>

\boldsymbol{a}(\boldsymbol{x}) = \frac{1}{\zeta } \overline{\boldsymbol{n}_{\mu }(\boldsymbol{x},\boldsymbol{y})}_{S_\mu }; \quad \quad \| \boldsymbol{a} \| =1; \quad \quad \zeta (\boldsymbol{x}) = \| \overline{\boldsymbol{n}_{\mu }(\boldsymbol{x},\boldsymbol{y})}_{S_\mu } \|     

 

</math>

where <math display="inline">\zeta (\boldsymbol{x})</math> is a measure of the ''tortuosity'' of the activated microscale failure path <math display="inline">S_{\mu }</math> (for instance, for <math display="inline">S_{\mu }</math> being a straight line, then <math display="inline">\zeta </math>=1). A similar expression for    this tortuosity factor can be found in Toro_al_FOMF_2014. In this case, the inelastic strain rate (Eq. [[#eq-2.36|2.36]]) is expressed as:

{\dot{\boldsymbol{\varepsilon }}}^{(i)}(\boldsymbol{x},t) = \zeta (\dot{\boldsymbol{\beta }}_{\mu }(\boldsymbol{x},t) \otimes \boldsymbol{a}(\boldsymbol{x}))^s     

 

</math>

In summary, Eq. [[#eq-2.35|2.35]] and Eq. [[#eq-2.36|2.36]] retrieve the format of a ''constitutive model equipped with an internal length and with internal variables whose evolution is  described by the microstructure behavior''. Although this model will never be used for computational purposes<span id="fnc-4"></span>[[#fn-4|¹]], it supplies relevant insights on the properties of the resulting homogenized  constitutive model.

<span id="fn-4"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-4|¹]]) Instead, the homogenized value of the stress in Eq. [[#eq-2.10|2.10]]  is point-wise used to evaluate the current macroscopic stress in terms of the corresponding macroscopic strain.</span>

In virtue of the Hill-Mandel Principle of Macro-Homogeneity, See:

 

 

Paper <math>\# </math>2''' it can be concluded that the macroscopic fracture energy <math display="inline">{G^{f}}</math> is equivalent to the average  of microscopic fracture energy <math display="inline">{G_{\mu }^{f}}(\boldsymbol{y})</math>, along the activated failure mechanism at the microscale <math display="inline">S_{\mu }</math>. Replacing Eq. [[#eq-2.41|2.41]] into Eq. [[#eq-2.9|2.9]], and after some  manipulations, the macroscopic fracture energy is given by the expression oliver2015continuum:

*  At the macroscale level, an EFEM based on the CSDA is used, aiming at capturing propagating cracks. As described in oliver2014crack, this technology consists of the insertion,  during specific stages of the simulation, of goal oriented specific strain fields via mixed finite element formulations. This allows modeling propagating cracks through the macroscale finite element  mesh.

* A standard Finite Element model is used at the microscale level, combining standard elements endowed with continuum hardening constitutive models and cohesive-band elements endowed with regularized constitutive softening models.  These are placed in the edges of every finite element, capturing the crack onset and strain localization, similar to the cohesive interface elements in pandolfi1999finite, and more  recently in rodrigues20162d. This approach benefits the simplicity of the algorithm and the non-intrusive character of its implementation.

|[[Image:draft_Samper_355780031-figure5.png|360px|]]

|[[Image:draft_Samper_355780031-figure6.png|210px|Multiscale model: finite element discretization at the microscale.]]

| style="text-align: center;" | <math>L_{\mu }=\sum _{\forall {\mathcal{B}^{(i)}_{\mu ,\hbox{act}}}\subset {\mathcal{B}_{\mu ,\hbox{coh}}}(\boldsymbol{y},t)} L_{\mu }^{(i)}; \quad \quad  {l_\mu }\equiv \frac{{\Omega _\mu }}{L_{\mu }} = \mathcal{O}({h_\mu }) </math>

At the microscale, where the morphology and the position of candidate propagation mechanisms are predefined, the two issues are of minor relevance due to the adopted  simplified failure-bands model.  However, at the macroscale, there is not a predefined failure path, and in principle, any material point may fail and propagate in any direction. To adequately solve  the previous questions, the procedure for modeling onset and propagation of discontinuities recently developed for monoscale problems  oliver2014crack has been extended to the multiscale setting. The proposed methodology is based on the use of the following specific techniques:

|[[Image:draft_Samper_355780031-figure10.png|480px|Evolution of the injection domains for three typical stages (t₁< t₂< t₃) of the  discontinuity propagation.]]

* '''Strain injection techniques''': based on the use of goal oriented assumed-strain fields injected in  selected domains, via mixed formulations Simo_Hughes_1986,Simo_Riffai_1990,Reddy_Simo_1995,zienk_taylor_2000. The standard (four points) Gauss quadrature rule, corresponding  to full integration of two-dimensional quadrilaterals, is complemented with two additional sampling points placed at the barycenter of the element (see Fig. [[#img-9|9]]),  termed ''singular'' and ''regular'' sampling points. These two additional quadrature points sample the stresses similarly to the standard Gauss points. Therefore, for the  injected elements, numerical integration (typically evaluation of the incremental internal forces in terms of the stresses), is based on those two additional sampling points by defining  the weight indicated in '''Table 3''' - in Appendix B in Paper <math display="inline">\#{2}</math>.

 

As for propagation purposes, two different  enhanced strain injection stages, are considered<span id="fnc-5"></span>[[#fn-5|¹]]:

 

<div id='img-9'></div>

|[[Image:draft_Samper_355780031-figure7.png|300px|Sampling points involved in the numerical integration.]]

\mathcal{D}(\boldsymbol{x},t) := \frac{1}{{\Omega _\mu }} \int _{{\mathcal{B}_{\mu }}} \mathcal{D}_\mu (\boldsymbol{y},t) \, dV = \int _{{\mathcal{B}_{\mu }}} ({\boldsymbol{\sigma }}_{\mu }: {\dot{\boldsymbol{\varepsilon }}_\mu }-\dot{\psi }_{\mu }) \, d{\mathcal{B}_{\mu }}\ge 0   

 

</math>

The onset of homogenized dissipation <math display="inline">\mathcal{D}(\boldsymbol{x},t)</math>, is used as the criterion to inject the weak discontinuity mode. This subset (''weak discontinuity domain'' <math display="inline">{\mathcal{B}_{\hbox{wd}}}</math>)<span id="fnc-6"></span>[[#fn-6|²]] is defined as:

{\mathcal{B}_{\hbox{wd}}}:=\left\{\bigcup _{e} {\mathcal{B}}^{(e)}; \quad {\mathcal{B}}^{(e)} \subset {\mathcal{B}_{\hbox{loc}}}(t); \quad 0<\mathcal{D}(\boldsymbol{x}_{C}^{(e)},t)<\mathcal{D}(\boldsymbol{x}_{C}^{(e)},t_{SD}) \right\}   

 

</math>

Eq. [[#eq-2.49|2.49]] characterizes <math display="inline">{\mathcal{B}_{\hbox{wd}}}</math> as the set of elements of the localization domain <math display="inline">{\mathcal{B}_{\hbox{loc}}}</math> whose barycenter has not yet bifurcated<span id="fnc-7"></span>[[#fn-7|³]],  and its corresponding homogenized dissipation take values different from zero, (and lower than <math display="inline">\mathcal{D}(\boldsymbol{x}_{C}^{(e)},t_{SD})</math>) in inelastic loading regime.

{\dot{\boldsymbol{\varepsilon }}}^{(e)}(\boldsymbol{x},t) \equiv \underbrace{\sum \limits _{i=1}^{n_{\hbox{node}}} {\boldsymbol{\nabla }}N_{i}(\boldsymbol{x}) \otimes {\dot{\boldsymbol{u}}}_{i}(t)}_{\hbox{regular}} + \underbrace{{\zeta _{\mathcal{S}}^{h^{(e)},l_{\mu }^{(e)}}}(\boldsymbol{x}){\dot{\boldsymbol{\gamma }}}^{(e)}(t)}_{\hbox{singular}}     

 

</math>

where <math display="inline">N_{i}</math> are the standard shape functions, <math display="inline">{\dot{\boldsymbol{u}}}(t)</math>, the macroscale nodal displacements, and <math display="inline">{\zeta _{\mathcal{S}}^{h^{(e)},l_{\mu }^{(e)}}}</math> is the regularized ''dipole-function'' in the element <math display="inline">(e)</math>. A  description of the weak enhanced mode is presented in Fig. [[#img-10|10]], In addition, the variational problem for the weak-discontinuity regime, in rate form, is presented in  '''Box A1''' in the Appendix B in Paper <math display="inline">\#{2}</math>.

 

<div id='img-10'></div>

|[[Image:draft_Samper_355780031-figure8.png|360px|Weak discontinuity mode. Elemental regularized dipole function ζ_\mathcalS^{h^(e),l_μ^(e)}.]]

* In a ''second stage, the strong discontinuity stage'', the obtained crack path field, <math display="inline">S</math>, is used to determine the appropriate placement of an elemental embedded strong  discontinuity strain field, which is incrementally injected in the appropriated set of elements <math display="inline">{\mathcal{B}_{\hbox{sd}}}</math>.

{\mathcal{B}_{\hbox{wd}}}:=\left\{\bigcup _{e} {\mathcal{B}}^{(e)}; \quad {\mathcal{B}}^{(e)} \subset {\mathcal{B}_{\hbox{loc}}}(t); \quad \mathcal{D}(\boldsymbol{x}_{C}^{(e)},t) \ge \mathcal{D}(\boldsymbol{x}_{C}^{(e)},t_{SD}) \right\}   

 

</math>

{\dot{\boldsymbol{u}}}= {\dot{\bar{\boldsymbol{u}}}}+ {\mathcal{H}_{\mathcal{S}}}\lbrack\lbrack{\dot{\boldsymbol{u}}}\rbrack\rbrack   

 

</math>

where <math display="inline">{\dot{\bar{\boldsymbol{u}}}}</math> stands for the smooth part of the displacement field, and <math display="inline">{\mathcal{H}_{\mathcal{S}}}</math> is the Heaviside function, shifted to <math display="inline">\mathcal{S}</math>. <p>

{\dot{\boldsymbol{\varepsilon }}}^{(e)}(\boldsymbol{x},t) \equiv \underbrace{\sum \limits _{i=1}^{\hbox{node}} {\boldsymbol{\nabla }}N_{i}(\boldsymbol{x}) \otimes {\dot{\boldsymbol{u}}}_{i}(t)}_{\hbox{regular}} + \underbrace{{\delta ^{l_{\mu }^{(e)}}_{s}}(\boldsymbol{\dot{\beta }}^{(e)} \otimes {\boldsymbol{n}}^{(e)})^{S}}_{\hbox{singular}}   

 

</math>

|[[Image:draft_Samper_355780031-figure9.png|300px|Strong discontinuity mode. Elemental regularized Dirac delta function δ^{l_μ^(e)}ₛ.]]

 

* '''Crack-path-field techniques''': their goal is the identification of the trace of the propagating crack by means of the so-called crack-path field.  It is denoted as <math display="inline">\mu (\boldsymbol{x},t)</math>, and obtained from a selected localized scalar variable <math display="inline">\alpha (\boldsymbol{x},t)</math>, which identifies the crack path as the locus where <math display="inline">\alpha (\boldsymbol{x},t)</math> takes its  transversal maximum  value. In order to define this locus <math display="inline">\mathcal{S}_t</math>, some alternatives have been developed in this multiscale framework, see  Eqs. <math display="inline">55-56</math> in Paper <math display="inline">\#{2}</math> oliver2015continuum. The variational  statement for the crack-path field model is fully detailed in '''Box 3.1''' in Paper <math display="inline">\#{2}</math>.

In regards to the space and time integrations, as commented above, injection of weak-discontinuity and strong-discontinuity modes requires, in principle, specific integration rules in space, : a standard four-point Gauss quadrature rule, and two additional sampling points, for injected elements, <math display="inline">{\mathcal{B}}^{(e)} \subset {\mathcal{B}_{\hbox{wd}}}</math> and <math display="inline">{\mathcal{B}}^{(e)} \subset {\mathcal{B}_{\hbox{sd}}}</math> so that <math display="inline">{\mathcal{B}_{\hbox{loc}}}= {\mathcal{B}_{\hbox{wd}}}\cup {\mathcal{B}_{\hbox{sd}}}</math>. Since those domains evolve along time (see Fig. [[#img-8|8]]), some additional problems on the time-integration of the  resulting equilibrium equations are found. To tackle this issue, in oliver2014crack and oliver2015continuum is proven that defining some "equivalent" stress entities  at the standard Gauss points, the spatial integration can be rephrased as a standard four Gauss points integration rule in the integration domain. This space-time integration rule is fully  explained in Appendix B3 in Paper <math display="inline">\#{2}</math>, and the corresponding stress evaluation is also summarized in Box A3.

<span id="fn-5"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-5|¹]]) To switch between stages, a set of ''control variables'' are defined, all those detailed in Sec. 3 and Appendix B in Paper <math>\#{2}</math></span>

<span id="fn-6"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-6|²]]) Under the  CSDA, the homogenized dissipation is evaluated at the barycenter of the finite element, denoted by <math>\boldsymbol{x}_{C}^{(e)}</math></span>

<span id="fn-7"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-7|³]]) For a deeper review of the bifurcation  analysis, and, the definition of the corresponding bifurcation time <math>t_{B}</math>, the reader is referred to Sec. 2.5 in Paper <math>\#{2}</math>, and, for numerical aspects Oliver_bifurc_2010.</span>

It can be immediately noticed the high computational cost demanded by the multiscale fracture model described in the previous chapter. This computational cost  becomes  unaffordable, even using resources of high-performance computing (also exhibited in the classical FE<math>^2</math> approach for smooth problems).

* Second Stage: also called ''HyPer-Reduced Order Model'' (HPROM). Two different techniques have been developed here. The first one is based on ''interpolation methods'',  widely applied in problems exhibiting hardening behavior (see Paper <math display="inline">\#{1}</math>) hernandez2014high. The second one, based on a Reduced  Optimal Quadrature (ROQ) rule, has been applied to fracture (non-smooth) problems. Similarly to the previous stage, the reduced modes functions used to reconstruct the  state variables, are computed off-line.

* Second Stage: is the development of the ''HyPer-Reduced Order Model'' (HPROM), and consist of redefining an adequate quadrature rule , demanding low computational cost. Two alternative quadrature rules have been developed. The first one is based on ''interpolation methods'', widely applied in problems exhibiting hardening behavior, and the second one, based on a Reduced  (''sub-optimal'') Quadrature rule (ROQ), being applied, in most cases, to fracture (non-smooth) problems. In the same way that the previous stage, the modes used to reconstruct the  state variables, are computed in the off-line part.

the latter takes advantage of the unbounded character of the strain field after the bifurcation time, originally intended for fracture processes, it can be  extended to smooth problems, being more general than the former one.

The model order reduction concept relies on the premise that, for any input parameter <math display="inline">\boldsymbol{\mu } \in \mathcal{D}</math> governing the microscale displacement fluctuations <math display="inline">{\boldsymbol{\tilde{u}}_\mu }</math>, the solution can be approximated by a set of <math display="inline">n</math> linearly independent basis functions <math display="inline">\boldsymbol{\Phi }</math> approximately spanning the ''primal variable''<span id="fnc-8"></span>[[#fn-8|¹]] space.

<span id="fn-8"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-8|¹]]) ''Primal variable'' is known as the selected  variable to perform the reduction process.</span>

Taking as a primal variable the displacement fluctuations, and departing from the problem depicted in Sec. [[#2.1 Computational Homogenization|2.1]], a first step consists of determining an approximation<span id="fnc-9"></span>[[#fn-9|¹]]  of the finite element space  of kinematically admissible microscale displacement fluctuations <math display="inline">\mathcal{\tilde V}^{h}_{\mu }</math>. This approximation is obtained as the span of the displacement fluctuation solutions obtained, for a judiciously  chosen set of <math display="inline">n_{hst}</math> input strain ''trajectories'', every trajectory being discretized into a number of steps <math display="inline">n_{stp}</math>. These set of finite element solutions are stored into the snapshot  matrix <math display="inline">\mathbf{X}_{u}</math> as column vectors:

Once the snapshot matrix <math display="inline">\mathbf{X}_{u}</math> has been computed, the ''Elastic-Inelastic decomposition'' technique is used to determine the ''reduced basis functions''.  The reason for it  relies on the fact that the SVD applied to the whole matrix <math display="inline">\mathbf{X}_{u}</math>, may produce basis with a large number of elements, which makes difficult to retrieve the response of the RVE in some specific cases. Particularly, the elastic response<span id="fnc-10"></span>[[#fn-10|²]], might request a much larger number of basis functions, this translating into a significant waste of computational cost.

 

 

Paper <math>\# </math>1'''obtaining the reduced basis as the combination (spatial sum) of both sub-bases. An orthonormal basis for <math display="inline">\mathcal{\tilde V}_{u,el}^{snp}</math> is determined  by taking a low number of elastic snapshots (at a minimum, <math display="inline">n_{snp}^e=3</math> for 2D problems, <math display="inline">n_{snp}^e=6</math> for 3D problems), and computing the corresponding orthonormal basis.

| style="text-align: center;" | <math>[\boldsymbol{\Phi }] = [ \underbrace{\boldsymbol{\Phi }_{1},\boldsymbol{\Phi }_{2},\boldsymbol{\Phi }_{3}}_{\hbox{elastic modes basis}} \quad    \underbrace{\boldsymbol{\Phi }_{4},\boldsymbol{\Phi }_{5},\boldsymbol{\Phi }_{6},\cdots ,\boldsymbol{\Phi }_{n_u}}_{\hbox{"essential" inelastic modes basis}} ] </math>

| style="text-align: center;" | <math>\mathcal{\tilde V}_{u}^{*}= \mathcal{\tilde V}_{u,el}^{snp}\oplus \mathcal{\tilde V}_{u,inel}^{snp,*}= \hbox{span} \left\{\underbrace{\boldsymbol{\Phi }_{1},\boldsymbol{\Phi }_{2},\boldsymbol{\Phi }_{3}}_{\hbox{elastic modes}} \quad    \underbrace{\boldsymbol{\Phi }_{4},\boldsymbol{\Phi }_{5},\boldsymbol{\Phi }_{6},\cdots ,\boldsymbol{\Phi }_{n_u}}_{\hbox{"essential" inelastic modes}} \right\}\subseteq \mathcal{\tilde V}^{h}_{\mu } </math>

Placing the <math display="inline">m_e</math> elastic modes in the first <math display="inline">m_e</math> positions, followed by the ''essential''<span id="fnc-11"></span>[[#fn-11|³]] inelastic modes, ensures the reduced-order model to deliver linear elastic solutions with the same accuracy than the HF solutions. For more details, the reader is encouraged to sent to the Appendix B  in hernandez2014high.

Once the reduced basis <math display="inline">[\boldsymbol{\Phi }]</math> is computed, See:

 

 

Paper <math>\# </math>1''' the online stage consists of solving the discrete version of the microscale  equilibrium equation (via FE), projected onto the reduced-order space <math display="inline">\mathcal{\tilde V}_{u}^{*}\subseteq \mathcal{\tilde V}^{h}_{\mu }</math> spanned by <math display="inline">[\boldsymbol{\Phi }]</math>. To this end, the test and trial functions, <math display="inline">{\boldsymbol{\eta }}</math> and <math display="inline">{\boldsymbol{\tilde{u}}_\mu }</math>, are approximated by the following linear expansions:

Introducing expressions [[#eq-3.6|3.6]] and [[#eq-3.7|3.7]] into the discrete version of the microscale BVP (see Sec. 4 in hernandez2014high), and multiplying the resulting expression by <math display="inline">\boldsymbol{\Phi }^{T}</math> (Galerkin projection), it yields:

minipage0.95 '''PROBLEM A (ROM)''' (''Microscale reduced problem via POD''):

Given the macroscale strain, }{\boldsymbol{\varepsilon }}\hbox{, and the reduced basis for displacement fluctuations }\boldsymbol{\Phi }\hbox{, find }\boldsymbol{c}\in \mathbb{R}^{n_{\varepsilon }}\hbox{ satisfying: equation  _B_ B^*T(y) _(y,+ B^*c,) dB__i=1^n_g B^*T(y_g,:) _(y_g,:) w_g = 0

where <math display="inline">\boldsymbol{c}=[c_1,c_2,\cdots ,c_{n_u}] \in \mathbb{R}^{n_u}</math> denotes the vector containing the coefficients associated to each basis function <math display="inline">\boldsymbol{\Phi }_i</math>, being <math display="inline">\boldsymbol{c}</math> the basic unknowns for the standard reduced-order  problem. <math display="inline">\mathbf{B}^*</math> stands for the ''reduced'' strain-displacement matrix “B-matrix” defined as <math display="inline">\mathbf{B}^{*}(\boldsymbol{y})=\mathbf{B}(\boldsymbol{y}) \cdot \boldsymbol{\Phi }(\boldsymbol{y})</math>. When using a Gauss quadrature  integration scheme, <math display="inline">n_g=\mathcal{O}(n)</math> is the total number of Gauss points of the mesh; <math display="inline">w_g</math> denotes the weight associated to the g-th Gauss point <math display="inline">\boldsymbol{y}_g</math>; <math display="inline">\mathbf{B}(\boldsymbol{y}_g,:)</math> and <math display="inline">{\boldsymbol{\sigma }}_{\mu }(\boldsymbol{y}_g,:)</math>  stand for the reduced B-matrix and the stress vector at Gauss point <math display="inline">\boldsymbol{y}_g</math>, respectively hdezMONOGR_141_HROM.

light-gray  '''Remark:''' The general procedure is performed into two parts, first, the ''offline part'' is devoted to compute the reduced basis functions in [[#3.2.1 Computation of the reduced basis functions|3.2.1]], including the sampling program to store the snapshots selected from pre-selected training trajectories, second, the ''online part'' is devoted to solve the PROBLEM A, see Eq. [[#eq-2.18|2.18]]. The  homogenization of the stress and constitutive tensors are also included in the ''online part''.

light-gray  '''Remark:''' However, the HPROM based on the interpolation of the stress field, has been proved to be effective in hardening processes, the Reduced Optimal Quadrature method (ROQ), is a more  general and robust strategy, it can be applied not only for hardening but also for softening processes. For details about the interpolation method, the reader is sent to read Secs. 5 to 7 in  hernandez2014high.

<span style="text-align: center; font-size: 75%;">([[#fnc-9|¹]]) In general, <math>\mathcal{\tilde V}^{h}_{\mu }</math> cannot be precisely determined, such a task will require  finite element analyses of the cell under ''all conceivable'' strain paths. Rather, one has to be content to construct an ''approximation'' of it.</span>

<span id="fn-10"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-10|²]]) Under an infinitesimal strain framework, this response is exactly recovered with only three basis hernandez2014high.</span>

<span id="fn-11"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-11|³]]) Essential based on a threshold given by an ''a-priori'' error estimation,  see Sec. 9.4 in hernandez2014high, thus, <math>\mathcal{\tilde V}_{u,inel}^{snp,*}</math> corresponds to the truncated version of the full base with <math>n_u-3</math> dominant modes.</span>

Taking advantage of the unbounded character of the microscale strain field typically observed in this kind of problems, the failure cell is splitted into a ''regular domain'' (made of elastic  matrix and possible inclusions) and a ''cohesive domain'' (cohesive bands exhibiting a softening cohesive behavior). Details on this issue can be found in Sec. 3.2.2 in Paper <math display="inline">\#{3}</math> Oliver_Caicedo_HROM_2017.

<span id="eq-3.8"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.8)

complemented by the evolution laws of the internal variables Simo_Hughes_1998. Thus, in consonance with the hierarchical multiscale approach, the following parametrized  functional can be defined:

<span id="eq-3.9"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.9)

minipage0.95  '''PROBLEM II (HF)''' (''Microscale saddle-point problem)'': Given the macroscale strain, }{\boldsymbol{\varepsilon }}\hbox{, find }{\boldsymbol{\tilde{\varepsilon }_{\mu }}}\hbox{ and }{\boldsymbol{\lambda }}\hbox{ satisfying: equation  _(,),(,) =arg\Big\{min__E_ _S^n n _[,] (_,)\Big\}  Such that: equation* =f(_,)

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.10)

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.11)

Eqs. [[#eq-3.10|3.10]] and [[#eq-3.11|3.11]] provide the solution of the saddle-point problem stated in Eq. [[#3.2.3 Formulation of the microscale saddle-point problem|3.2.3]]. It can be proven that Eqs.  [[#eq-3.10|3.10]] and [[#eq-3.11|3.11]] make ''PROBLEM II'' equivalent to the original problem in Eq. [[#eq-2.11|2.11]], but now rephrased in terms of the microscale strain fluctuations <math display="inline">{\boldsymbol{\tilde{\varepsilon }_{\mu }}}</math> (see ''PROBLEM I-R'' in Paper <math display="inline">\#{3}</math>).

<span id="eq-3.12"></span>

<span id="eq-3.12"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.12)

where <math display="inline">n_{snp}</math> is the number of snapshots vectors. See:

 

 

Paper <math>\# </math>3''' Therefore, <math display="inline">\mathbf{X}</math> represents a number of sampled solutions obtained with the HF model  under different loading conditions. For more details, the reader is addressed to Sec. 4 in Oliver_HROM_2017.

|[[Image:draft_Samper_355780031-figure4.png|540px|(a) All entries of \mathbfX are partitions into two sub-blocks: \mathbfX_reg and \mathbfX_coh, the strain fluctuations from points    at ''regular'' and ''cohesive'' domains, respectively. (b) Snapshots taken from the elastic regime of the failure cell correspond to the sub-block \mathbfX^E. Snapshots taken during the    inelastic regime (at least one Gauss point at \mathcalB_μ,\hboxcoh is in inelastic state) correspond to the sub-block \mathbfX^I.]]

In addition, the ''Elastic-Inelastic snapshot decomposition'' above explained hernandez2014high, See:

'''Sec.3.2.2

 

 

<span id="eq-3.13"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.13)

where <math display="inline">\mathbf{X}_{\hbox{reg}}</math> and <math display="inline">\mathbf{X}_{\hbox{coh}}</math> collect the strain fluctuations located outside and inside the cohesive bands, respectively. The right hand side matrix  in [[#eq-3.13|3.13]] emphasizes the double partition performed in accordance with elastic-inelastic regimes.

With [[#eq-3.13|3.13]], the next step consists on computing the reduced basis corresponding to the elastic behavior of each subdomain. Those basis are computed via SVD, and  this technique is applied, in a separate way, for each partition corresponding to <math display="inline">{\mathcal{B}_{\mu }}\setminus {\mathcal{B}_{\mu ,\hbox{coh}}}</math> and <math display="inline">{\mathcal{B}_{\mu ,\hbox{coh}}}</math>, resulting in the following reduced basis:

<span id="eq-3.14"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.14)

<span id="eq-3.15"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.15)

<span id="eq-3.16"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.16)

where each element <math display="inline">{\boldsymbol{\Psi }}_i</math>, of the basis <math display="inline">{\boldsymbol{\Psi }}_G</math>, is a microscale strain fluctuation mode and the vector of time dependent coefficients <math display="inline">\boldsymbol{c}(t)=[c_1,c_2,\cdots ,c_{n_{\boldsymbol{\varepsilon }}}]</math> (<math display="inline">\boldsymbol{c} \in \mathbb{R}^{n_{\boldsymbol{\varepsilon }}}</math>) represents their corresponding amplitudes (the actual unknowns of the problem). See:  

 

 

<span id="eq-3.17"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.17)

minipage0.95  '''PROBLEM III (ROM)''' (RVE ''saddle point problem)'':  Given the macroscale strain, }{\boldsymbol{\varepsilon }}\hbox{, find }\mathbf{c}\in \mathbb{R}^{n_{\varepsilon }}\hbox{ and }{\boldsymbol{\lambda }}\in \mathbb{R}^{n_{\sigma }}\hbox{ satisfying: align c(,),(,)  = arg\Big\{min_cR_n_ _R^n_ _[,] ( c,)\Big\}  equation* with _[,] ( c,)= _B_ _(+ c) dB_+ ^T \left(_B_  dB_\right)c  such that:  equation* =f(_,)

<span id="eq-3.18"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.18)

<span id="eq-3.19"></span>

<span id="eq-3.20"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.19)

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.20)

<span id="eq-3.21"></span>

<span id="eq-3.22"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.21)

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.22)

where <math display="inline">[{\boldsymbol{\sigma }}_{\mu }(\boldsymbol{c})]_{G} \in \mathbb{R}^{(n_{{\boldsymbol{\sigma }}} \cdot N_{g})}</math> is the column vector constituted by piling-up the <math display="inline">N_{g}</math> stress vectors, <math display="inline">{\boldsymbol{\sigma }}_{\mu }(\boldsymbol{c}) \in \mathbb{R}^{(n_{{\boldsymbol{\sigma }}})}</math>, evaluated at the integration Gauss points. The column vector <math display="inline">[{\boldsymbol{\lambda }}]_{G}</math> is also the pilled-up of <math display="inline">N_g</math> repeated values of the same constants vector <math display="inline">{\boldsymbol{\lambda }}\in \mathbb{R}^{(n_{{\boldsymbol{\sigma }}})}</math>.  See:  

 

 

<span id="eq-3.23"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.23)

light-gray  Considering the system of equations [[#eq-3.21|3.21]] and [[#eq-3.22|3.22]] for <math>\boldsymbol{c}</math> and <math>{\boldsymbol{\lambda }}</math>, it could be expected that this problem, of <math>n_{{\boldsymbol{\varepsilon }}} + n_{\sigma }</math> equations, should be less computationally costly, than the HF model. However, this is not the case: the numerical simulations with the ROM model do not  substantially reduce the computational cost, and little (or even smaller than one) speedups are obtained. This fact highlights that the actual bottleneck for fast online computation is not the  solution of the balance equations but, rather, the determination of the stresses, internal forces and stiffness matrices at every integration point of the underlying finite element mesh. Therefore, an additional technique is proposed to reduce the amount of integration points in which the constitutive equation is evaluated.

light-gray  '''Remark (B):''' The general procedure is performed into two parts, first, the ''offline part'' is devoted to compute the reduced basis functions in [[#3.2.3.1 Computation of the reduced basis functions|3.2.3.1]], including the  sampling program to store the snapshots selected from pre-selected training trajectories, second, the ''online part'' is devoted to solve the PROBLEM III, see Eq. [[#3.2.3.2 Formulation of the reduced order model|3.2.3.2]].  The homogenization of the stress and constitutive tensors are also included in the ''online part''.

Next, an additional reduction step to diminish the computational burden for computing equation [[#eq-2.18|2.18]]  or [[#eq-3.21|3.21]] and [[#eq-3.22|3.22]] is introduced.  In addition, in order to guarantee the good performance of the numerical integration procedure, all possible ''a-priori'' operators<span id="fnc-12"></span>[[#fn-12|¹]]  have to be computed during the offline part. Particularly, the term <math display="inline">[W]\cdot [{\boldsymbol{\Psi }}_{G}]</math> in Eq.  [[#eq-3.22|3.22]] can be computed entirely in the offline part.

<span id="fn-12"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-12|¹]]) ''A-priori'' operators can be interpreted as all  possible algebraic operations that remain constant during the computation and that can be gathered in a ''global'' operator, computed during the offline stage. These operators are stored and called during the online stage for the global assembling  of the discrete version of the optimality conditions (equilibrium equations).</span>

The problem of constructing a <math display="inline">\mathcal{O}(n_u)</math>-dimensional representation of the stress field is similar to the procedure previously presented for reducing the dimensionality of the displacement  fluctuation field <math display="inline">{\boldsymbol{\tilde{u}}_\mu }</math>. The goal is to find a set of orthogonal basis functions, such that its span accurately approximates the set of all possible ''statically admissible stress solutions''.  In accordance with the procedure described in Sec. sec:SmoothPr_ROM, first, the finite element stress distribution for representative macroscale strain histories are computed. Then,  the ''Elastic-Inelastic decomposition'' (see. Sec. 3.2.1 in hernandez2014high) is applied to the resulting set of stress solutions <math display="inline">[{\boldsymbol{\sigma }}_{\mu }^1(\boldsymbol{y}),{\boldsymbol{\sigma }}_{\mu }^2(\boldsymbol{y}),\dots ,{\boldsymbol{\sigma }}_{\mu }^{n_{stp}}(\boldsymbol{y})]</math> in order to identify both the elastic and the ''essential'' inelastic stress modes.

| style="text-align: center;" | <math>\mathcal{V}_{\sigma }^{*}= \hbox{span}\left\{\underbrace{{\boldsymbol{\Psi }}_1(\boldsymbol{y}),{\boldsymbol{\Psi }}_2(\boldsymbol{y}),{\boldsymbol{\Psi }}_3(\boldsymbol{y}),\dots ,{\boldsymbol{\Psi }}_{m_{e}}(\boldsymbol{y})}_{\hbox{Elastic stress modes}},\underbrace{{\boldsymbol{\Psi }}_{m_{e}+1},{\boldsymbol{\Psi }}_{m_{e}+2}\dots ,{\boldsymbol{\Psi }}_{n_{{\boldsymbol{\sigma }}}}}_{\hbox{"Essential" inelastic stress modes}}\right\} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.24)

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.25)

<span id="eq-3.26"></span>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3.26)

'''Remark:''' A close examination of Eq. [[#eq-3.26|3.26]] reveals that, the "standard" strategy proves completely fruitless, for it leads to ''ill-posed''  reduced order model equations. The reason relies on the fact that, every element of the stress field snapshot matrix <math display="inline">\mathbf{X}_{\boldsymbol{\sigma }}</math> is, by construction, obtained by a lineal combination of elements  contained in <math display="inline">{\boldsymbol{\Psi }}</math>,  thus, in consonance with Eq. [[#eq-2.18|2.18]], in which the representative cell does not experience external forces,  the space spanned by the reduced basis <math display="inline">{\boldsymbol{\Psi }}</math> is contained into the null space of the operator <math display="inline">\mathbf{B}^{*T}</math>, namely, See:

'''Sec. 5.2