# A priori hyper-reduction method for coupled viscoelastic-viscoplastic composites

Bilel Miled, David Ryckelynck, Sabine Cantournet

► To cite this version:

Bilel Miled, David Ryckelynck, Sabine Cantournet. A priori hyper-reduction method for coupled viscoelastic-viscoplastic composites. Computers & Structures, 2013, 119, pp.95-103. 10.1016/j.compstruc.2012.11.017. hal-00790175

HAL Id: hal-00790175

https://minesparis-psl.hal.science/hal-00790175

Submitted on 16 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A priori hyper-reduction method for coupled viscoelastic–viscoplastic composites

B. Miled, D. Ryckelynck, S. Cantournet

Centre des matériaux, UMR CNRS 7633, Ecole des mines de Paris, France

In this paper, a mean field homogenization (MFH) method is compared to the hyper reduction (HR) method. The homogenization of concern aims to forecast the mechanical response of viscoelastic visco-plastic composites undergoing small strains. Reference results are provided by the usual finite element method (FEM) applied to an unit cell problem. In both methods the microscopic strain fields are represented using a reduced basis. In MFH it is an eigenstrain basis in the vocabulary of [17]. This basis is spanned by the stress free strains introduced by Eshelby [5]. In the HR method the reduced basis is spanned by modes. It can be created by the proper orthogonal decomposition (POD) method or the APHR method [19]. MFH and HR methods are compared in terms of equation formulation, accuracy and computational time. The accuracy of both global and local results are compared. We consider as MFH local results the global ones, as if they are uniform in the matrix of the composite. It turns out that the HR method provides simulations of accuracy and computational complexity between the MFH method and the full field FEM. The HR model contains a reduced mesh named reduced domain (RD). This requires to reconstruct the internal variables by using the Gappy POD. We point out that the APHR method provides unrealistic non-smooth modes when the reconstruction of the internal variables is performed only outside the RD and not inside the RD.

**Keywords:** Reduced order model, Reduced domain Finite elements, Viscoelasticity, Viscoplasticity
# 1. Introduction

Polymer matrix composites reinforced with glass fibers are widely used in a variety of technological applications because of the enhanced mechanical properties, particularly the stiffness and the strength. Among these applications, one can cite the air craft wings, some car components (e.g. body frames, hood and door panels), and the sporting equipment. The principal difficulties associated with predicting the mechanical behavior of these com posites are, on the one hand, the complexity of the geometry (e.g. volume fraction, shape and orientation of glass fibers) and on the other hand the constitutive behavior of the thermoplastic polymer matrix (e.g. rate dependent, hardening, temperature). Predicting of the effective response of viscoelastic viscoplastic thermoplastic polymers reinforced with linear elastic spherical particles requir ing low computational cost and good accuracy is the main illustra tion of this paper.

Accurate predictions of the macroscopic mechanical response of such composite materials may be derived from full field calcula tion (e.g. by finite element methods) of the local stresses and strains throughout a statistically representative volume element (RVE) of the microstructure subjected to periodic boundary condi tions. Such modeling is however too costly when applied to the simulation of real scale structures made of composite materials. Alternatively, various multiscale computational strategies have been developed to tackle this important issue, such as homogeni zation strategies (e.g. mean field homogenization (MFH), asymp totic or mathematical theory of homogenization, methods of cells and subcells associated with the transformation fields analysis), or reduced order models (ROMs) (a posteriori and a priori meth ods) which can efficiently be relied upon. Their remarkably low computational cost as compared to full field computation over a RVE allows them to be used as implicit constitutive laws in large scale simulations.

The present contribution is concerned with the application of ROM and MFH methods to coupled viscoelastic viscoplastic com posites. In both methods the microscopic strain fields are repre sented using a reduced basis (RB). In MFH it is an eigenstrain basis in the vocabulary of [17]. This RB is spanned by the stress free strains introduced by Eshelby [5]. In the HR method the re duced basis is spanned by modes. It can be created by the POD method or the APHR method [19]. MFH and HR methods are com pared in terms of equation formulation, accuracy and computa tional time. The accuracy of both global and local results are compared. We consider as MFH local results the global ones, as if they are uniform in the matrix of the composite.

In the scientific literature, many authors shown that the recourse to ROMs can speed up the solution of complex finite element simu lations. The aim of these models is to supply few shape functions to represent the spatially distributed state of a given system. In non linear time dependent problems, the most used methods to extract a reduced basis from a set of fields provided by simulations using classical approximation methods (i.e. finite element method) are the POD and the snapshot POD. The following papers give access to relevant literature for both approaches: [12,13,23,24] and [25]. The POD method was proposed to solve homogenization problem in [26] and [16]. But using a POD Galerkin formulation of the equilibrium equation do not reduce the computational complexity of the internal variable computation. The details about this draw back are presented in Section 2.2. An alternative, approach to POD method is the non incremental proper generalized decomposition (PGD) method. This method was proposed to solve homogenization problem in [9,10] and [3]. Non incremental stands for iterative time integration schemes that provide prediction over the entire time interval at each iteration. As explained in [22], this kind of

method is not appropriate to all constitutive equation in mechanics of materials. Therefore, an incremental approach based on usual time integration schemes has been developed. This approach is the hyper reduction method [19] extended to mechanical model involving internal variables in [20]. Assuming that constitutive equations are local equations, the HR method introduces a reduced domain to state the weak form of the equilibrium equations related to a given RB. By reducing the integration domain, some explanatory internal variables are estimated through the constitutive equation related to the RD. These explanatory variables are the internal vari ables of the RD. The remainder of the internal variables, outside the RD, is explained by the explanatory variables by solving least square minimization problem as proposed in the Gappy POD [6]. This pro cedure requires a RB dedicated to the internal variables. In this pa per, all the RB are provided by the APHR method. This method is an adaptive approach that enrich the RBs during the approximate solu tion of the FE equations.

Concerning the homogenization strategy applied in this work namely MFH method, we note that this technique is based on the interaction laws between the different phases. It provides an approximative effective response, as well as a description of mechanical fields within the phases in terms of volume averages, requiring a low computational cost. The mean field methods were first proposed for composites having linear elastic or thermo elastic constituents. The generalization of this homogenization procedure for nonlinear two phase composites is based on the lin earization of the local constitutive laws and the per phase unifor mization of the mechanical properties. Those two conditions being satisfied, the nonlinear homogenization problem becomes similar to an elastic or a thermo elastic homogenization problem, depend ing on the choice of the linearization process. This kind of problem enable to introduced stress free strains as proposed by Eshelby [5]. For the viscoelastic viscoplastic behavior, only the so called incre mentally affine formulation was applied. This method has been proposed for elasto viscoplastic composites [4], and after general ized for coupled VE VP behavior [15]. This proposal is valid for multi axial, non monotonic and non proportional loading histo ries, and leads to thermoelastic like relations directly in the time domain. This linearization formulation is also valid for any visco plastic model.

The paper is organized in the following manner: Section 2 intro duces the coupled viscoelastic viscoplastic constitutive law, the formulation of the equilibrium equations related to the HR method and the APHR algorithm used to generate the RB. The mean field homogenization method based on an incrementally affine lineari zation method is summarized in Section 3. The comparison be tween the numerical results obtained with the FE method, HR method, and the MFH method are discussed in Section 4. Section 5 is the conclusion of this paper.

The following abbreviations are used throughout the text. VE: viscoelastic(ity), VP: viscoplastic(ity), EVP: elasto viscoplastic(ity), VE VP: viscoelastic(ity) viscoplastic(ity), MFH: mean field homog enization, RB: reduced basis, ROM: reduced order models, APHR: a priori hyper reduction, POD: proper orthogonal decomposition, and FE: finite elements.

Boldface symbols designate second or fourth rank tensors, as indicated by the context. Dyadic and inner products are expressed as:

$$ (\mathbf{a} \otimes \mathbf{b})_{ijkl} \quad a_{ij}b_{kl}, \quad \mathbf{a} : \mathbf{b} \quad a_{ij}b_{ji}, \quad (\mathbf{A} : \mathbf{b})_{ij} \quad A_{ijkl}b_{lk} $$

where summation over a repeated index is supposed. The symbols **1** and **I** designate the second and fourth rank symmetric identity tensors, respectively. Finally, the spherical and deviatoric operators **I**<sup>vol</sup> and **I**<sup>dev</sup> are given by:

$$ \mathbf{I}^{\text{vol}} \equiv \frac{1}{3} \mathbf{1} \otimes \mathbf{1} \text{ and } \mathbf{I}^{\text{dev}} \equiv \mathbf{I} \quad \mathbf{I}^{\text{vol}} $$
so that for $a_{ij} = a_{ji}$ we have:

$$ \mathbf{I}^{\text{vol}} : \mathbf{a} = \frac{1}{3} a_{mm} \mathbf{1}, \quad \mathbf{I}^{\text{dev}} : \mathbf{a} = \mathbf{a} = \frac{1}{3} a_{mm} \mathbf{1} \equiv \text{dev}(\mathbf{a}) $$

## 2. Formulation of the hyper-reduced model

### 2.1. The reference model

In this work, a coupled VE VP constitutive law proposed by Miled et al. [14] is considered. For this model, the total strain is supposed to be the sum of VE and VP parts, and the Cauchy stress is related to the history of VE strains via a linear VE model written as a Boltzmann integral [2]. The general form of the considered problem is given by this system, for all point x in the domain $\Omega$:

$$ \left\{ \begin{array}{l} \text{div}\boldsymbol{\sigma} \quad \mathbf{0} \\ \boldsymbol{\epsilon} \quad \boldsymbol{\epsilon}^{\text{ve}} + \boldsymbol{\epsilon}^{\text{vp}} \\ \boldsymbol{\sigma}(t) \quad \int_{\infty}^{t} \mathbf{C}^{\text{ve}}(t \ \tau) : \frac{\partial \boldsymbol{\epsilon}^{\text{ve}}}{\partial \tau} d\tau \\ \dot{\boldsymbol{\epsilon}}^{\text{vp}} \quad \dot{p} \frac{\partial f}{\partial \boldsymbol{\sigma}} \\ \dot{p} \quad g_{\nu}(\boldsymbol{\sigma}_{\text{eq}}, p, \dot{\boldsymbol{\epsilon}}) \ge 0 \end{array} \right. \quad (1) $$

For an isotropic material, the fourth order relaxation tensor is written as:

$$ \mathbf{C}^{\text{ve}}(t) = 2G(t)\mathbf{I}^{\text{dev}} + 3K(t)\mathbf{I}^{\text{vol}} \quad (2) $$

where $G(t)$ and $K(t)$ are shear and bulk relaxation functions, respectively, which can be expressed in the form of Prony series:

$$ \begin{aligned} G(t) = & G_{\infty} + \sum_{i=1}^{l} G_i \exp \left( \frac{t}{g_i} \right) \quad \text{and} \quad K(t) \\ & + \sum_{j=1}^{J} K_j \exp \left( \frac{t}{k_j} \right). \end{aligned} \quad (3) $$

Here, $g_i$ ($i = 1...I$) and $k_j$ ($j = 1...J$) are the deviatoric and volumetric relaxation times respectively; $G_i$ ($i = 1...I$) and $K_j$ ($j = 1...J$) are the corresponding moduli or weights, and $G_{\infty}$ and $K_{\infty}$ are the long term elastic shear and bulk moduli.

The yield function $f$ is expressed in function of the von Mises equivalent stress $\sigma_{eq}$, the initial yield stress $\sigma_y$ (which may depend on the strain rate) and the hardening stress $R(p)$ (which depend on the accumulated plastic strain $p$):

$$ f(\sigma_{eq}, p, \dot{\boldsymbol{\epsilon}}) = \sigma_{eq} \quad (\sigma_y(\dot{\boldsymbol{\epsilon}}) + R(p)); \quad \sigma_{eq} \equiv \sqrt{\frac{3}{2}} \mathbf{s} : \mathbf{s} \quad (4) $$

For this work, we consider a power law hardening function:

$$ R(p) = k p^n \quad (5) $$

where $k$ [MPa] is the hardening modulus and $n$ the hardening exponent. A power law VP function is also considered with two parameters: the viscoplastic modulus ($\kappa[1/s]$) and exponent ($m$) which appear as follows:

$$ g_{\nu}(\sigma_{\text{eq}}, p) = \begin{cases} \kappa \left( \frac{f}{\sigma_{y}+R(p)} \right)^{m} & \text{if } f > 0 \\ 0 & \text{otherwise} \end{cases} \quad (6) $$

The deviatoric and dilatational parts of the stress tensor can be expressed in function of the deviatoric and dilatational parts of the strain tensor using Eqs. (1) (3):

$$ \left\{ \begin{array}{l} \mathbf{s}(t) = 2G_{\infty} \xi^{\text{ve}}(t) + \sum_{i=1}^{l} \underbrace{2G_i \exp \left( \frac{t}{g_i} \right) \int_{\infty}^{t} \exp \left( \frac{\tau}{g_i} \right) \frac{\partial \xi^{\text{ve}}}{\partial \tau} d\tau}_{\mathbf{s}_i(t)} \\ \sigma_H(t) = 3K_{\infty} \epsilon_H^{\text{ve}}(t) + \sum_{j=1}^{J} \underbrace{3K_j \exp \left( \frac{t}{k_j} \right) \int_{\infty}^{t} \exp \left( \frac{\tau}{k_j} \right) \frac{\partial \epsilon_H^{\text{ve}}}{\partial \tau} d\tau}_{\sigma_{H_j}(t)} \end{array} \right. \quad (7) $$

In order to compute the previous integral equations, two time integration methods have been tested. The first one supposes that the VE strain rate is constant in a generic time interval $[t_n, t_{n+1}]$. The second method is the mid point integration rule. A return mapping algorithm has been also proposed based on two steps, VE predictor followed by VP corrector. After some algebra, the following total stress computed at $t_{n+1}$ has been found:

$$ \boldsymbol{\sigma}(t_{n+1}) = \boldsymbol{\sigma}^{\text{pred}}(t_{n+1}) = \tilde{\mathbf{C}}^{\text{ve}} : \Delta\boldsymbol{\epsilon}^{\text{vp}} \quad (8) $$

with

$$ \begin{aligned} \boldsymbol{\sigma}^{\text{pred}}(t_{n+1}) = & \mathbf{C}_{\infty}^{\text{ve}} : \boldsymbol{\epsilon}^{\text{ve}}(t_n) + \tilde{\mathbf{C}}^{\text{ve}} \\ & : \Delta\boldsymbol{\epsilon} + \sum_{i=1}^{l} \exp \left( \frac{\Delta t}{g_i} \right) \mathbf{s}_i(t_n) \\ & + \sum_{j=1}^{J} \exp \left( \frac{\Delta t}{k_j} \right) \boldsymbol{\sigma}_{H_j}(t_n) \mathbf{1} \end{aligned} \quad (9) $$

where $\mathbf{C}_{\infty}^{\text{ve}} = 2G_{\infty}\mathbf{I}^{\text{dev}} + 3K_{\infty}\mathbf{I}^{\text{vol}}$ and $\tilde{\mathbf{C}}^{\text{ve}} = 2\tilde{\mathbf{G}}^{\text{dev}} + 3\tilde{\mathbf{K}}^{\text{vol}}$. The increase in mental relaxation moduli $\tilde{\mathbf{G}}$ and $\tilde{\mathbf{K}}$ which are functions of the time increment, are expressed for the first integration method as follows:

$$ \left\{ \begin{array}{l} \tilde{\mathbf{G}} \equiv G_{\infty} + \sum_{i=1}^{l} G_i \left[ 1 \exp \left( \frac{\Delta t}{g_i} \right) \frac{g_i}{\Delta t} \right] \\ \tilde{\mathbf{K}} \equiv K_{\infty} + \sum_{j=1}^{J} K_j \left[ 1 \exp \left( \frac{\Delta t}{k_j} \right) \frac{k_j}{\Delta t} \right] \end{array} \right. \quad (10) $$

and for the second integration method:

$$ \tilde{\mathbf{G}} \quad G\left(\frac{\Delta t}{2}\right) \quad \text{and} \quad \tilde{\mathbf{K}} \quad K\left(\frac{\Delta t}{2}\right) \quad (11) $$

### 2.2. The hyper reduced approximation

#### 2.2.1. Weak form based on a reduced domain

Usually, in mechanical engineering, model reduction methods aim to introduce a RB to represent the displacement. Let us denote $(\phi_k^u)_{k=1...M_u}^u$ the set of modes of the reduced basis. One have to fore cast the linear combination of the vectors of the RB according to a weak form of Problem (1). The usual Galerkin procedure uses the vectors of the RB as test functions to state the following weak form of the mechanical problem: find $(\phi_k^u(t))_{k=1}^{M_u}$, $\epsilon^{\text{ve}}$, $\epsilon^{\text{vp}}$, $p$ and $\boldsymbol{\sigma}$ such that:

$$ \left\{ \begin{array}{l} \mathbf{u}(\mathbf{x}, t) = \sum_{k=1}^{M_u} \phi_k^u(\mathbf{x}) \phi_k^u(t) \\ \int_{\Omega} \phi_k^u \text{div}\boldsymbol{\sigma} \, d\Omega = \mathbf{0} \quad \forall k = 1, \dots, M_u \\ \boldsymbol{\epsilon}(\mathbf{u}) = \boldsymbol{\epsilon}^{\text{ve}} + \boldsymbol{\epsilon}^{\text{vp}} \quad \forall \mathbf{x} \in \Omega \\ \boldsymbol{\sigma}(t) = \int_{\infty}^{t} \mathbf{C}^{\text{ve}}(t \ \tau) : \frac{\partial \boldsymbol{\epsilon}^{\text{ve}}}{\partial \tau} d\tau \quad \forall \mathbf{x} \in \Omega \\ \dot{\boldsymbol{\epsilon}}^{\text{vp}} = \dot{p} \frac{\partial f}{\partial \boldsymbol{\sigma}} \quad \forall \mathbf{x} \in \Omega \\ \dot{p} = g_{\nu}(\sigma_{\text{eq}}, p, \dot{\boldsymbol{\epsilon}}) \ge 0 \quad \forall \mathbf{x} \in \Omega \end{array} \right. \quad (12) $$

Using a numerical one step time integration scheme, one can fore cast different states of the system at different instants. According to an incremental formulation, for a generic time interval $[t_i, t_{i+1}]$, the mechanical state of the variables (stress field and the internal variables) is assumed to be known at time $t_i$. The unknowns are the state variables at time $t_{i+1}$. Obviously, the integral over the domain $\Omega$ requires the prediction of the stress field and the internal variables over the entire domain $\Omega$. This does not reduce the computational complexity related to the constitutive equations.

The HR method aims to introduce a RD, denoted $\Omega_Z$, by using truncated test functions $(\psi_k^u)_{k=1}^{M_u}$ in the Sobolev space function $H^1(\Omega)$. A truncated vector is associated to each vector of the RB re
lated to displacements. Inside the RD, each truncated mode $\psi_k^u$ matches the empirical mode $\phi_k^u$ such that the following con strained is fulfilled outside the RD:

$$ \psi_k^u(\mathbf{x}) \quad 0 \quad \forall \mathbf{x} \in \Omega \setminus \Omega_Z \quad \forall k \quad 1, \dots, M_u \qquad (13) $$

Therefore, we obtain the following weak form of the mechanical problem: find $(\phi_k^u(t))_{k, 1}^{M_u}, \hat{\epsilon}^{ve}, \hat{\epsilon}^{vp}, \hat{p}$ and $\hat{\sigma}$ such that:

$$ \left\{ \begin{array}{l} \mathbf{u} \sum_{k=1}^{M_u} \phi_k^u \phi_k^u(t) \\ \int_{\Omega_Z} \psi_k^u \mathrm{div} \hat{\sigma} \, d\Omega \quad \mathbf{0} \quad \forall k \quad 1, \dots, M_u \\ \epsilon(\mathbf{u}) \quad \hat{\epsilon}^{ve} + \hat{\epsilon}^{vp} \quad \forall \mathbf{x} \in \Omega_Z \\ \hat{\sigma}(t) \quad \int_{-\infty}^{t} \mathbf{C}^{ve}(t) \, \tau \, : \frac{\partial}{\partial t} \hat{\epsilon}^{ve} \, d\tau \quad \forall \mathbf{x} \in \Omega_Z \\ \hat{\epsilon}^{vp} \quad \hat{p} \frac{\partial \hat{p}}{\partial \hat{\sigma}} \quad \forall \mathbf{x} \in \Omega_Z \\ \hat{p} \quad g_v(\hat{\sigma}_{eq}, \hat{p}, \hat{\epsilon}) \ge 0 \quad \forall \mathbf{x} \in \Omega_Z \end{array} \right. \qquad (14) $$

Here the symbol $\hat{\cdot}$ is used for the variables restricted to $\Omega_Z$. These variables are named explanatory variables. The smaller $\Omega_Z$, the less complex are the set of constitutive equations in problem (14). But, $\Omega_Z$ must be large enough to have a well posed problem with a unique solution.

The mechanical problem (14) is well posed if the tangent stiff ness matrix rank is $M_u$. Here, the tangent stiffness matrix reads:

$$ K_{kj} \int_{\Omega_Z} \psi_k^u \mathrm{div}(\mathbf{C}^{\mathrm{alg}} : \epsilon(\phi_j^u)) \, d\Omega \qquad (15) $$

The fourth order tensor $\mathbf{C}^{\mathrm{alg}}$ is the algorithmic tangent operator (Eq. (29)) related to the local constitutive equation. If during the numerical simulation of (14) the tangent stiffness matrix is rank deficient, then the numerical simulation must be restarted with a bigger $\Omega_Z$. For instance, one can add the neighbour element of $\Omega_Z$ in the extended $\Omega_Z$. This procedure converge toward a Galerkin formulation when $\Omega_Z$ tends to $\Omega$. In practice a convenient extent of $\Omega_Z$ is directly provided by the method presented in Section 2.3.

### 2.2.2. Full field reconstruction

In Eq. (14) the internal variables and the stress are forecast only inside $\Omega_Z$. At this stage, it is not a full field prediction because these fields are not estimated outside $\Omega_Z$. We use the Gappy POD method to estimate the internal variables over $\Omega \setminus \Omega_Z$ by using a RB for each internal variable. For instance, when the RB dedicated to displace ments is created, we simultaneously generate a RB for accumulated plastic strains. This RB is denoted $(\phi_k^p)_{k, 1, M_p}$. The same procedure holds for all these variables: $\epsilon^{ve}, \epsilon^{vp}$ and $\sigma$. The Gappy POD reads: find $(\phi_k^p)_{k, 1, M_p}$ such that,

$$ \left\{ \begin{array}{l} p(\mathbf{x}, t) \quad \hat{p} \quad \forall \mathbf{x} \in \Omega_Z \\ p(\mathbf{x}, t) \quad \sum_{k=1}^{M_p} \phi_k^p(\mathbf{x}) \phi_k^p(t) \quad \forall \mathbf{x} \in \Omega \setminus \Omega_Z \\ (\phi_k^p)_{k, 1 \dots M_p} \quad \mathrm{argmin}_{(y_k)_{k=1}^{M_p}} \int_{\Omega_Z} \hat{p} \sum_{k=1}^{M_p} \phi_k^p y_k \quad d\Omega \end{array} \right. \qquad (16) $$

In the above extrapolation, $p$ is a repaired value of the accumulated plastic strain by using the Gappy POD. Here, we preserve inside $\Omega_Z$ the predictions provided by the constitutive equations. According to [6], it corresponds to the available value $\hat{p}$ at all points of $\Omega_Z$. As we will show in the following illustrations, unphysical non smooth field may be generated by this procedure. For instance, it is the case when the RB involves only one constant field and when $\hat{p}$ is not uni form over $\Omega_Z$. Therefore, we proposed a smooth Gappy POD. This extrapolation method reads:

$$ \left\{ \begin{array}{l} p(\mathbf{x}, t) \quad \sum_{k=1}^{M_p} \phi_k^p(\mathbf{x}) \phi_k^p(t) \quad \forall \mathbf{x} \in \Omega \\ (\phi_k^p)_{k, 1 \dots M_p} \quad \mathrm{argmin}_{(y_k)_{k=1}^{M_p}} \int_{\Omega_Z} \hat{p} \sum_{k=1}^{M_p} \phi_k^p y_k \quad d\Omega \end{array} \right. \qquad (17) $$

In (17) the field of $p$ is fully represented by the reduced basis $(\phi_k^p)_{k, 1}^{M_p}$. This prediction erases the one given by the local constitutive equation.

### 2.3. Determination of the $\Omega_Z$

The APHR method is an incremental method. At each time increment there is a predictor step and if the accuracy of the prediction is not satisfactory a corrector step is performed. The corrector step is a full FE solution, over one time increment, which is initialized with the results of the predictor step. It provides the corrections denoted $\delta\mathbf{u}, \delta\mathbf{p}, \delta\epsilon^{ve}, \delta\epsilon^{vp}, \delta\sigma$. These corrections are defined over $\Omega$. They are used as new vectors to expand the RB related to each of these variables. More details about this adaptive procedure can be found in [21].

The predictor step is performed by the HR method. But, if the prediction is not smooth despite the FE solution is smooth, then the correction is going to be non smooth. Therefore, the RB will in volve unphysical non smooth vectors. Hence, we prescribe to use the smooth Gappy POD to avoid the generation of unphysical RB.

At each time step $[t_i, t_{i+1}]$, four information are provided by the APHR method: (1) a prediction of the stress field and the internal variables over the entire domain $\Omega$; (2) convenient reduced bases $(\phi_k^u$ for displacement field and $\phi_k^p$ for accumulated plastic strain) to represent the prediction of the state variables over the time interval; (3) a reduced domain $\Omega_Z$; (4) truncated modes $\psi_k^u$ related to the reduced basis $\phi_k^u$.

The main idea of the hyper reduction is to write both equilibrium equations and local constitutive equations at points of $\Omega$ where the strains induced by the modes are large, and where the magnitude of the internal variables is large. A RB being known, the construction of $\Omega_Z$ follows five steps:

$$ \Omega_k^p \quad \mathrm{arg} \max_{\epsilon_k^{\Omega}, i, 1, \dots, N_e} \int_{\Omega_i^{\Omega}} (\phi_k^p)^2 \, d\Omega, \quad k, 1, \dots, M_p \qquad (18) $$

$$ \Omega_k^u \quad \mathrm{arg} \max_{\epsilon_k^{\Omega}, i, 1, \dots, N_e} \int_{\Omega_i^{\Omega}} \epsilon(\phi_k^u) : \epsilon(\phi_k^u) \, d\Omega, \quad k, 1, \dots, M_u \qquad (19) $$

$$ \bar{\Omega}^u \quad \cup_{k=1}^{M_u} \Omega_k^u \qquad (20) $$

$$ \tilde{\mathcal{F}}_Z \quad \{i \mid i \in \{1, \dots, N\}, \quad \|\mathbf{N}_i\|_{\bar{\Omega}^u} > 0\} \qquad (21) $$

$$ \Omega_Z \quad \{ \epsilon_j^\Omega \mid j \in 1, \dots, N_e, \quad \sum_{i \in \tilde{\mathcal{F}}_Z} \| \mathbf{N}_i \|_{\epsilon_j^\Omega} > 0 \} \quad \cup \quad \cup_{k=1}^{M_p} \Omega_k^p \qquad (22) $$

Here, $N$ is the number of FE shape functions, $\mathbf{N}_i$ is the $i$th FE shape function, $N_e$ is the number of elements of the mesh, and $\epsilon_j^\Omega$ is the support of the $j$th element. Eq. (18) enables to locate the elements of the mesh where the internal variables magnitude is large. Eq. (19) enables to locate the elements of the mesh where the strain magnitude is large. Eq. (21) gives the set of nodes connected to the set of elements covering the domain $\bar{\Omega}^u$, and Eq. (22) gives the union of the supports of the FE shape functions of the set of nodes $\tilde{\mathcal{F}}_Z$ completed by the contributions of the modes related to the internal variables. One can notice in Eqs. (20) and (22) that the bigger the RB, the bigger the RD. In practice, the extent and the shape of $\Omega_Z$ is sufficient to ensure that Eqs. (14), (16) and (17) have unique solutions. As mention above, The extent of the RD can be increased by adding to the list, the elements connected to at least one element of the list.
## 2.4. Error indicator

An error estimator is defined by using the partial residual de noted $\mathbf{R}_Z$ of the parametrized problem according to the Riesz rep resentation theorem. The definition of the full residual reads:

$$ \mathbf{R}^N \in \mathbb{R}/R_i \quad \int_{\Omega} N_i(x) \mathrm{div} \sigma \, d\Omega \quad \forall i \ 1, \dots, N \qquad (23) $$

Unfortunately, due to the hyper reduction, the constitutive equa tions provide a stress prediction only in $\Omega_Z$. Therefore, only a part of the equilibrium residuals are available. They are related to the degrees of freedom that are not connected to the remainder part $\Omega \setminus \Omega_Z$ of the domain. We denote $\mathcal{F}_Z$ the list of dof indexes that are not connected to $\Omega \setminus \Omega_Z$. Therefore the partial residual reads:

$$ \mathbf{R}_Z^N \in \mathbb{R}/R_{Z_i} \quad \int_{\Omega_Z} N_i(x) \mathrm{div} \hat{\sigma} \, d\Omega \quad \forall i \in \mathcal{F}_Z \qquad (24) $$

$$ R_{Z_i} = 0 \quad \text{if } i \notin \mathcal{F}_Z \qquad (25) $$

The definition of the error estimator, denoted by $\eta$, reads:

$$ \eta = \|\mathbf{R}_Z\| \qquad (26) $$

Obviously, this error indicator may under estimate the true error. An approximate prediction is satisfactory if:

$$ \eta \le \epsilon_R \qquad (27) $$

where $\epsilon_R$ is the reduced order model error.

# 3. Mean-field homogenization method

The mean field homogenization procedure for nonlinear two phase composites is based on the linearization of the local consti tutive laws and the per phase uniformization of the mechanical properties. The linearization formulation adopted in [15] for VE VP behavior with a coupling between linear isotropic VE and $J_2$ VP, under isothermal conditions and a small perturbation hypoth esis, is the incrementally affine linearization method. The basic equations and the main hypotheses of this linearization method are summarized in this section.

## 3.1. The incrementally affine formulation

In order to find the incrementally affine relation, we start by the linearization of the evolution equations of the viscoplastic strain and the accumulated plastic strain at the beginning of a time step around the end time of the step. Next, a numerical integration of the linearized equations is required using a fully implicit backward Euler scheme. The obtained algebraic equations lead to the follow ing incrementally affine formulation (to more detail, see [15]):

$$ \Delta\sigma \quad \mathbf{C}^{\mathrm{alg}}(t_{n+1}) : (\Delta\epsilon \quad \Delta\epsilon^{\mathrm{af}}) \qquad (28) $$

The algorithmic tangent operator at time $t_{n+1}$ ($\mathbf{C}^{\mathrm{alg}}(t_{n+1})$) is ex pressed as follows:

$$ \mathbf{C}^{\mathrm{alg}} \quad \tilde{\mathbf{C}}^{\mathrm{ve}} \quad \frac{(2\tilde{G})^2}{h_v} \mathbf{N} \otimes \mathbf{N} \quad (2\tilde{G})^2 \quad \frac{\sigma_{eq} \Delta p}{\sigma_{eq} + 3\tilde{G} \Delta p} \frac{\partial \mathbf{N}}{\partial \sigma} \quad \frac{2\tilde{G}}{h_v g_{,\sigma}} \mathbf{N} \otimes g_{,\epsilon} \qquad (29) $$

where the rate dependent denominator $h_v$ is defined by this expression:

$$ h_v = \frac{1}{(\Delta t) g_{,\sigma}} + 3 \tilde{G} \frac{g_{,p}}{g_{,\sigma}} \qquad (30) $$

$$ \text{Notations: } g_{,\sigma} = \frac{\partial g_{,v}}{\partial \sigma_{eq}}; \quad g_{,p} = \frac{\partial g_{,v}}{\partial p}; \quad g_{,\epsilon} = \frac{1}{\Delta t} \frac{\partial g_{,v}}{\partial \epsilon} $$

$\mathbf{N}$ is the normal to the yield surface in stress space, and the expres sion of the affine strain increment $(\Delta\epsilon^{\mathrm{af}})$ is:

$$ \Delta\epsilon^{\mathrm{af}} \quad \Delta\epsilon^{\mathrm{af}}_{\mathrm{evp}} + \frac{1}{2\tilde{G}} \sum_{i} \left[ 1 \exp\left(\frac{\Delta t}{g_{,i}}\right) \mathbf{s}_{i}(t_{n}) + \frac{1}{3\tilde{K}} \sum_{j} \left[ 1 \exp\left(\frac{\Delta t}{k_{j}}\right) \sigma_{H_{j}}(t_{n}) \mathbf{1} \right] \right] \qquad (31) $$

In $J_2$ viscoplasticity, the expression of the affine strain increment for the EVP behavior $\Delta\epsilon^{\mathrm{af}}_{\mathrm{evp}}$ is:

$$ \Delta\epsilon_{\mathrm{evp}}^{\mathrm{af}} \quad \dot{p}(t_n) \Delta t \left[ \mathbf{N}(t_n) + \mathbf{N}(t_{n+1}) \frac{g_p(t_{n+1}) \Delta t}{g_p(t_{n+1}) \Delta t} \right] \qquad (32) $$

## 3.2. Analogy with linear thermoelasticity

Starting from the incrementally affine expression (28), Miled [15] shown that this equation is form similar to linear thermoelastic relation ($\sigma = \mathbf{C}^{\mathrm{el}}; \epsilon + \beta\theta$), so that at each time step, available homogenization models for linear thermoelastic composites can be applied for nonlinear VE VP composites using the following substitutions:

$$ \sigma \rightarrow \Delta\sigma, \quad \mathbf{C}^{\mathrm{el}} \rightarrow \mathbf{C}^{\mathrm{alg}}, \quad \epsilon \rightarrow \Delta\epsilon, \quad (\beta\theta) \rightarrow \mathbf{C}^{\mathrm{alg}} : \Delta\epsilon^{\mathrm{af}} \qquad (33) $$

where $\theta$ designates a change in temperature and $\mathbf{C}^{\mathrm{el}}$ is the elastic stiffness modulus. Knowing the homogenization procedure of linear thermoelastic composites (see e.g. [11] and [18]), and using the sub stitutions of Eq. (33), the macroscopic stress tensor has been found as follows:

$$ \langle \Delta\sigma \rangle \quad \bar{\mathbf{C}} : \langle \Delta\epsilon \rangle \quad v_0 \hat{\mathbf{C}}_0^{\mathrm{alg}} : \Delta\hat{\epsilon}_0^{\mathrm{af}} \quad v_1 \hat{\mathbf{C}}_1^{\mathrm{alg}} : \Delta\hat{\epsilon}_1^{\mathrm{af}} \quad v_1 (\hat{\mathbf{C}}_1^{\mathrm{alg}} \quad \hat{\mathbf{C}}_0^{\mathrm{alg}}) : (\mathbf{A} \quad \mathbf{I}) \quad (\hat{\mathbf{C}}_1^{\mathrm{alg}} : \Delta\hat{\epsilon}_0^{\mathrm{af}} \quad \hat{\mathbf{C}}_0^{\mathrm{alg}} : \Delta\hat{\epsilon}_0^{\mathrm{af}}) \qquad (34) $$

where $\bar{\mathbf{C}}$ is the effective stiffness of the composite:

$$ \bar{\mathbf{C}} \quad [v_0 \hat{\mathbf{C}}_0^{\mathrm{alg}} + v_1 \hat{\mathbf{C}}_1^{\mathrm{alg}} : \mathbf{B}] : [v_0 \mathbf{I} + v_1 \mathbf{B}] \quad 1 \qquad (35) $$

$\mathbf{A}$ and $\mathbf{B}$ are the strain concentration tensors which are related as follows:

$$ \mathbf{A} \quad \mathbf{B} : (v_1 \mathbf{B} + v_0 \mathbf{I}) \quad 1 \qquad (36) $$

For this work, the strain concentration tensor of the Mori Tanaka scheme is used:

$$ \mathbf{B} \quad [\mathbf{I} + \mathbf{P} : (\hat{\mathbf{C}}_1^{\mathrm{alg}} \quad \hat{\mathbf{C}}_0^{\mathrm{alg}})] \quad 1 \qquad (37) $$

where $\mathbf{P} \quad \mathbf{S} : (\hat{\mathbf{C}}_0^{\mathrm{alg}})^{-1}$ is Hill's polarization tensor and $\mathbf{S}$ is Eshelby's tensor, which depends only on the properties of the matrix and the inclusion shape. We note that the comparison algorithmic tangent operators for the matrix $(\hat{\mathbf{C}}_0^{\mathrm{alg}})$ and inclusion $(\hat{\mathbf{C}}_1^{\mathrm{alg}})$ phases are uni form and anisotropic, by construction.

# 4. Numerical analysis

In order to assess the numerical simulations of both MFH and APHR methods, we have carried out FE simulations of a continuous VE VP matrix reinforced by spherical inclusions with linear elastic behavior. For this, the numerical algorithm introduced for the VE VP model was integrated into the ZeBuLon FE software [1] as a User defined MATerial (UMAT) routine.

The volume fraction of particles is 15%. Assuming a periodic microstructure and uniaxial stress tests, a 2D axisymmetric unit cell is defined and depicted in Fig. 1a. The prescribed boundary conditions are the following. Zero displacements are imposed in the radial direction on the left vertical side and in the vertical direction on the bottom horizontal side. The right vertical side is
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** 2D unit cell representation of a two-phase composite with periodic microstructure under axisymmetric loading: (a) full domain, (b) reduced domain in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 1**
Constitutive model parameters for polycarbonate at 22 °C.</caption><thead><tr><th colspan="5">Viscoelastic parameters</th></tr><tr><th>Initial shear modulus</th><th colspan="3">G<sub>0</sub> = 1074 MPa</th></tr><tr><th>Initial bulk modulus</th><th colspan="3">K<sub>0</sub> = 3222 MPa</th></tr><tr><th>G<sub>i</sub> (MPa)</th><th>g<sub>i</sub> (s)</th><th>K<sub>i</sub> (MPa)</th><th>k<sub>i</sub> (s)</th></tr></thead><tbody><tr><td>157</td><td>0.0021</td><td>472</td><td>0.007</td></tr><tr><td>80</td><td>0.00378</td><td>242</td><td>0.126</td></tr><tr><td>37</td><td>0.0248</td><td>111</td><td>0.216</td></tr><tr><th colspan="5">Viscoplastic parameters</th></tr><tr><th>Hardening function</th><td>k = 150 MPa</td><td colspan="3">n = 0.43</td></tr><tr><th>Viscoplastic function</th><td>κ = 150 s<sup>-1</sup></td><td colspan="3">m = 5</td></tr><tr><th>Yield stress</th><td>σ<sub>y</sub> = 35 MPa</td><td colspan="3"></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Finite-element, incrementally affine homogenization, and a priori hyper reduction results of the stress–strain curves in tension test at the strain rate of 1 s⁻¹ for a volume fraction of inclusions equal to 15%.</figcaption></figure>

<table><caption>**Table 2**
Accuracy of each computational method computed at the end of the simulation.</caption><thead><tr><th></th><th>j = MFH</th><th>j = REF1</th><th>j = REF2</th><th>j = REF3</th><th>j = REF4</th></tr></thead><tbody><tr><th>|d<sub>i</sub>/d<sub>FE</sub>| (%)</th><td>8.5</td><td>0.03</td><td>0.14</td><td>0.67</td><td>0.67</td></tr><tr><th>|p<sub>i</sub>/p<sub>FE</sub>| (%)</th><td>5.22</td><td>0.08</td><td>5.72</td><td>0.21</td><td>0.21</td></tr><tr><th>|P<sub>max</sub>/P<sub>max</sub>| (%)</th><td>79.64</td><td>0.1</td><td>0.86</td><td>0.1</td><td>0.09</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Average accumulated plastic strain in the matrix after tensile deformation up to 5% for a volume fraction of the inclusions equal to 15% and at the strain rate of 1 s⁻¹.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Maximum accumulated plastic strain in the matrix after tensile deformation up to 5% for a volume fraction of the inclusions equal to 15% and at the strain rate of 1 s⁻¹.</figcaption></figure>
<table><caption>**Table 3**
Computational times.</caption><thead><tr><th></th><th>FE</th><th>MFH</th><th>REF1</th><th>REF2</th><th>REF3</th><th>REF4</th></tr></thead><tbody><tr><td>CPU time (s)</td><td>407</td><td>4</td><td>62</td><td>74</td><td>46</td><td>156</td></tr><tr><td>CPU time saving (%)</td><td>0</td><td>99</td><td>85</td><td>81</td><td>89</td><td>61</td></tr></tbody></table>

constrained to have uniform radial displacement. Uniform vertical displacement is imposed on the top horizontal side. We used the GMSH [8] software to mesh the geometry, and a typical mesh comprises approximately 914 elements and 5198 number of degrees of freedom (dof). The FE computations are conducted using six noded triangular elements (CAX6) for the inclusion phase and eight noded quadratic elements(CAX8) for the matrix phase.

The elastic properties of the particles are $E_1 = 76$ GPa and $v_1 = 0.22$, and those of the matrix are collected in Table 1 which presents parameters of polycarbonate at 22 °C based on experimental data collected from [7].

In this section, different strategies are adopted to compute the numerical results. The first strategy consists in the FE simulation

which represents the reference solution. The second strategy is related to the MFH method. The third and fourth strategies introduce the adaptive procedure with hyper reduction method using two different approaches to extrapolate the field of internal variables: a smooth Gappy POD (17) for the third strategy and the usual Gappy POD (16) for the fourth strategy. The fifth strategy is reproducing the third without adaptation of the ROM. It is a pure HR simulation based on the smooth Gappy POD. The sixth strategy is reproducing the third without adaptation of the ROM and without hyper reduction. It is the usual REF4alerkin strategy. These six strategies are denoted respectively FE strategy (FE), semi analytical strategy (MFH), APHR strategies (REF1, REF2, REF3), and POD Galerkin strategy (REF4). We note that for the strategies REF1 and REF2, the reduced state variables are computed without using an initial set of known shape functions. However, the RB obtained with the REF1 strategy is used to compute the reduced state variables for the strategies REF3 and REF4.

The RD consists of all elements in red that can be seen in Fig. 1b. Its mesh comprises approximately 133 elements.

The curves stress strain obtained for the reference simulation (FE), reduced model (APHR), and the MFH method are shown in

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** modes $(\epsilon_{22}(\psi_k^H))_{k=1}^{M_u}$ found with: (a) REF1 strategy and (b) REF2 strategy.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Modes of accumulated plastic strain field found with: (a) REF1 strategy and (b) REF2 strategy.</figcaption></figure>

Fig. 2. Compared to the FE analysis, the macroscopic predictions of the APHR method are very good whatever the adopted strategy (the accuracy of each computational strategy is reported in Table 2). However, the correlation between FE and mean field simulations is rather low (the discrepancy between both predictions is of the order of 8.5%).

In order to investigate the accuracy of the predictions of the APHR method to capture the plastic strains, Figs. 3 and 4 show the development of average and maximum accumulated plastic strains, respectively, in the matrix surrounding the spherical inclusion. Each figure includes the results of three different APHR strategies. The conclusion that can be drawn, in terms of correlation between APHR and FE methods, is that the REF3 strategy gives very good results for all cases. However, the results of the REF1 strategy are good for the average accumulated plastic strain, and the results of the REF2 strategy are good for the maximum accumulated plastic strain.

The computational time related to each strategy is reported in Table 3. For the APHR method, the CPU time saving varies between 61% for the REF4 strategy and 89% for the REF3 strategy. This gain can be contributed at the same time to the reduction of the size of the problem (i.e. the number of vectors of the reduced basis), and

to the hyper reduction, because thanks to the hyper reduction, the calculation at each time step is not performed on all elements of the integration domain, but only on a restrained number of elements. This comparison between the computational times shows also that the MFH method has a bigger contribution to the saving of the computational time compared to the APHR method.

The observation of modes which constitute the reduced basis (used for the APHR computation) allows a better understanding of the considered problem. In practice, in a POD reduced basis, the first modes are representative of global modes. Then the higher modes are associated with localized phenomena. This remark is illustrated in Fig. 5 where we present the different modes corresponding to the second component of the deformation field obtained from the APHR model using the REF1 strategy (six modes) and the REF2 strategy (seven modes). In Fig. 6, we can find the different modes corresponding to the accumulated plastic strain field obtained from the APHR model using the REF1 strategy (five modes) and the REF2 strategy (eight modes).

It is important to note that the number of modes forming the reduced basis depends strictly on the chosen HR strategy. Here, we can conclude that the non smooth HR strategy provides non smooth
and un physical modes both for strains and accumulated plastic
strain. To preserve the same order of accuracy, the non smooth HR
strategy requires more modes than the smooth HR strategy.

5. Conclusions

In this work, we studied a technique for the construction of re-
duced model from a finite element model for coupled VE VP solids.
It is based on the use of the APHR method which was previously
introduced by Ryckelynck [19] in the framework of thermome-
chanical simulations and has been generalized here for VE VP
composites involving material non linearities. In opposition with
the a posteriori decomposition techniques, the APHR method is a
priori resolution technique which is driving to a POD without
requiring any initial ROM nor known solutions. This technique is
based on two attractive features namely the adaptive and the hyper
reduction strategies. For the first feature, the reduced basis repre
senting all the previous numerical results obtained during the sim
ulation is updated during each new increment. For the second one,
the recourse to the hyper reduction of the equations provides sig
nificant computational time savings without losing much in terms
of result accuracy. In this work, various strategies are tested by
using usual Gappy POD or smooth Gappy POD during the hyper
reduction. In comparison with the solution obtained for the full
model with classical FE methods, the computation time necessary
for solving the numerical example treated is reduced, whereas the
accuracy is of the same order of magnitude. However, in compari
son with the predictions of the MFH method, the computation time
is more expensive, whereas the accuracy is improved.

The obtained results are overall encouraging for the application
of the APHR method in the coupled VE VP case. However, modeling
of the mechanical behavior of composite materials in the nonlinear
domain merits further developments. First, the performances of the
APHR method should be validated using more advanced nonlinear
models, in particular VE VP behavior with nonlinear VE part or
nonlinear kinematic hardening VP part. Another research direction
is to study the solution of large scale parametric problem by using a
multidimensional formulation. Finally, The study of nonlinear
behavior could be accompanied by a study of damage phenomena
(inclusion fracture, interface decohesion, matrix damage, etc.)

Acknowledgements

This work is related to the CSDL and CALME projects.

References

[1] Besson J, Foerch R. Large scale object-oriented finite element code design. Comput Methods Appl Mech Eng 1997;142:165-87.

[2] Boltzmann L. Wied Ann 1878;5:430.

[3] Chinesta F, Ladeveze P, Cueto E. A short review on model order reduction based on proper generalized decomposition. Arch Comput Methods Eng 2011;18(4):395-404.

[4] Doghri I, Adam L, Bilger N. Mean-field homogenization of elasto-viscoplastic composites based on a general incrementally affine linearization method. Int J Plasticity 2010;26:219-38.

[5] Eshelby JD. The determination of the elastic field of an ellipsoid inclusion, and related problems. Proc Roy Soc London A 1957;241:376-96.

[6] Everson R, Sirovich L. Karhunen Loève procedure for gappy data. J Opt Soc Am A 1995;12:1657-64.

[7] Frank GJ. Analytic and experimental evaluation of the effects of temperature and strain rate on the mechanical response of polymers. Report UDR-TR-97-152 1997, University of Dayton Research Institute, Dayton, Ohio.

[8] Geuzaine C, Remacle JF. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int J Numer Methods Eng 2009;79:1309-31.

[9] Ladevèze P, Passieux JC, Neron D. The LATIN multiscale computational method and the proper generalized decomposition. Comput Methods Appl Mech Eng 2010;199(21-22):1287-96.

[10] Lamari H, Ammar A, Cartraud P, Legrain G, Chinesta F, Jacquemin F. Routes for efficient computational homogenization of nonlinear materials using the proper generalized decompositions. Arch Comput Methods Eng 2010;17(4):373-91.

[11] Lielens G. Micro-macro modeling of structured materials. PhD thesis, Université catholique de Louvain, Belgium; 1999.

[12] Lorenz EN. Empirical orthogonal functions and statistical weather prediction. Scientific report 1, MIT Department of Meteorology, Statistical Forecasting Project; 1956.

[13] Lumley J. The structure of inhomogeneous turbulence. Atmos Turbul Wave Propag 1967;166-78.

[14] Miled B, Doghri I, Delannay L. Coupled viscoelastic-viscoplastic modeling of homogeneous and isotropic polymers: numerical algorithm and analytical solutions. Comput Methods Appl Mech Eng 2011;200:3381-94.

[15] Miled B. Coupled viscoelastic-viscoplastic modeling of homogeneous and reinforced thermoplastic polymers. PhD thesis, Université catholique de Louvain, Belgium; 2011.

[16] Monteiro E, Yvonnet J, He QC. Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction. Comput Mater Sci 2008;42(4):704-12.

[17] Mura T. Micromechanics of defects in solids. Dordrecht: Martinus Nijhoff Publishers; 1987.

[18] Pierard O, Friebel C, Doghri I. Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation. Compos Sci Technol 2004;64:1587-603.

[19] Ryckelynck D. A priori hyperreduction method: an adaptive approach. J Comput Phys 2005;202:346-66.

[20] Ryckelynck D. Hyper reduction of mechanical models involving internal variables. Int J Numer Methods Eng 2009;77(1):75-89.

[21] Ryckelynck D, Missoum Benziane D. Multi-level a priori hyper-reduction of mechanical models involving internal variables. Comput Methods Appl Mech Eng 2010;199:1134-42.

[22] Ryckelynck D, Vincent F, Cantournet S. Multidimensional a priori hyper-reduction of mechanical models involving internal variables. Comput Methods Appl Mech Eng 2012;225-228:28-43.

[23] Sirovich L. Quart Appl Math XLV 1987;3:561-71.

[24] Sirovich L. Turbulence and the dynamics of coherent structures Part I: Coherent structures. Quart Appl Math 1987;65:561-71.

[25] Sirovich L. Turbulence and the dynamics of coherent structures. Part III: Dynamics and scaling. Quart Appl Math 1987;65:583-90.

[26] Yvonnet J, Zahrouni H, Potier-Ferry M. A model reduction method for the post-buckling analysis of cellular microstructuresfla. Comput Methods Appl Mech Eng 2007;197(104):265-80.