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

In this paper, a mean ﬁeld 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 ﬁnite element method (FEM) applied to an unit cell problem. In both methods the microscopic strain ﬁelds are repre sented 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 com putational 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 ﬁeld 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 pro vides unrealistic non smooth modes when the reconstruction of the internal variables is performed only outside the RD and not inside the RD.


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.
Boldface symbols designate second or fourth rank tensors, as indicated by the context.Dyadic and inner products are expressed as: ða bÞ ijkl a ij b kl ; a : b a ij b ji ; ðA : bÞ ij 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 vol and I dev are given by:  [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 X: For an isotropic material, the fourth order relaxation tensor is writ ten as: where G(t) and K(t) are shear and bulk relaxation functions, respec tively, which can be expressed in the form of Prony series: Here, g i (i = 1 . . .I) and k j (j = 1 . . .J) are the deviatoric and volumet ric relaxation times respectively; G i (i = 1 . . .I) and K j (j = 1 . . .J) are the corresponding moduli or weights, and G 1 and K 1 are the long term elastic shear and bulk moduli.
The yield function f is expressed in function of the von Mises equivalent stress r eq , the initial yield stress r y (which may depend on the strain rate) and the hardening stress R(p) (which depend on the accumulated plastic strain p): For this work, we consider a power law hardening function: where k [MPa] is the hardening modulus and n the hardening expo nent.A power law VP function is also considered with two param eters: the viscoplastic modulus (j[1/s]) and exponent (m) which appear as follows: The deviatoric and dilatational parts of the stress tensor can be ex pressed in function of the deviatoric and dilatational parts of the strain tensor using Eqs.(1) (3): In order to compute the previous integral equations, two time inte gration 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 sec ond 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: where The incre mental relaxation moduli e G and e K which are functions of the time increment, are expressed for the first integration method as follows: and for the second integration method: uðx; tÞ 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 vari ables) 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 X requires the prediction of the stress field and the internal vari ables over the entire domain X.This does not reduce the computa tional complexity related to the constitutive equations.
The HR method aims to introduce a RD, denoted X Z , by using truncated test functions w u k À Á Mu k 1 in the Sobolev space function H 1 (X).A truncated vector is associated to each vector of the RB re lated to displacements.Inside the RD, each truncated mode w u k matches the empirical mode / u k such that the following con strained is fulfilled outside the RD: Therefore, we obtain the following weak form of the mechanical problem: find u u k ðtÞ À Á Mu k 1 ; ve ; vp ; p and r such that: Here the symbol ^is used for the variables restricted to X Z .These variables are named explanatory variables.The smaller X Z , the less complex are the set of constitutive equations in problem (14).But, X 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: The fourth order tensor C alg is the algorithmic tangent operator (Eq.( 29)) related to the local constitutive equation.If during the numer ical simulation of ( 14) the tangent stiffness matrix is rank deficient, then the numerical simulation must be restarted with a bigger X Z .For instance, one can add the neighbour element of X Z in the ex tended X Z .This procedure converge toward a Galerkin formulation when X Z tends to X.In practice a convenient extent of X Z is directly provided by the method presented in Section 2.3.

Full field reconstruction
In Eq. ( 14) the internal variables and the stress are forecast only inside X Z .At this stage, it is not a full field prediction because these fields are not estimated outside X Z .We use the Gappy POD method to estimate the internal variables over XnX 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 accumu lated plastic strains.This RB is denoted / p k À Á k 1...Mp .The same procedure holds for all these variables: ve , vp and r.The Gappy POD reads: find In the above extrapolation, p is a repaired value of the accumulated plastic strain by using the Gappy POD.Here, we preserve inside X Z the predictions provided by the constitutive equations.According to [6], it corresponds to the available value p at all points of X 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 p is not uni form over X Z .Therefore, we proposed a smooth Gappy POD.This extrapolation method reads: pðx; tÞ In (17) the field of p is fully represented by the reduced basis This prediction erases the one given by the local constitu tive equation.

Determination of the X Z
The APHR method is an incremental method.At each time increment there is a predictor step and if the accuracy of the pre diction is not satisfactory a corrector step is performed.The correc tor step is a full FE solution, over one time increment, which is initialized with the results of the predictor step.It provide the cor rections denoted du, dp, d ve , d vp , dr.These corrections are de fined over X.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 X; (2) convenient reduced bases (/ u k for displacement field and / p k for accumulated plastic strain) to represent the prediction of the state variables over the time inter val; (3) a reduced domain X Z ; (4) truncated modes w u k related to the reduced basis / u k .The main idea of the hyper reduction is to write both equilib rium equations and local constitutive equations at points of X 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 X Z follows five steps: Here, N is the number of FE shape functions, N i is the ith FE shape function, N e is the number of elements of the mesh, and e X j is the support of the jth 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 con nected to the set of elements covering the domain X u , and Eq. ( 22) gives the union of the supports of the FE shape functions of the set of nodes e 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 ex tent and the shape of X 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 con nected to at least one element of the list.

Error indicator
An error estimator is defined by using the partial residual de noted R Z of the parametrized problem according to the Riesz rep resentation theorem.The definition of the full residual reads: R N 2 R=R i Z X N i ðxÞ divr dX 8i 1; . . .; N ð 23Þ Unfortunately, due to the hyper reduction, the constitutive equa tions provide a stress prediction only in X 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 XnX Z of the domain.We denote F Z the list of dof indexes that are not connected to XnX Z .Therefore the partial residual reads: The definition of the error estimator, denoted by g, reads: Obviously, this error indicator may under estimate the true error.An approximate prediction is satisfactory if: where R is the reduced order model error.

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.

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]): The algorithmic tangent operator at time t n+1 (C alg (t n+1 )) is ex pressed as follows: where the rate dependent denominator h v is defined by this expression: Notations : g ;r @g v @r eq ; g ;p @g v @p ; g ; 1 Dt N is the normal to the yield surface in stress space, and the expres sion of the affine strain increment (D af ) is: In J 2 viscoplasticity, the expression of the affine strain increment for the EVP behavior D af evp is:

Analogy with linear thermoelasticity
Starting from the incrementally affine expression (28), Miled [15] shown that this equation is form similar to linear thermoelas tic relation (r = C el : + bh), 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: where h designates a change in temperature and C 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: where C is the effective stiffness of the composite: A and B are the strain concentration tensors which are related as follows: For this work, the strain concentration tensor of the Mori Tanaka scheme is used: where P S : b C alg 0 1 is Hill's polarization tensor and 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 b C alg 0 and inclusion b C alg 1 phases are uni form and anisotropic, by construction.

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   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 1 for a volume fraction of inclusions equal to 15%.

Table 2
Accuracy of each computational method computed at the end of the simulation.j = MFH j = REF1 j = REF2 j = REF3 j = REF4  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 com prises 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 m 1 = 0.22, and those of the matrix are collected in Table 1 which presents parameters of polycarbonate at 22 °C based on experi mental 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 re lated 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 reduc tion.It is the usual REF4alerkin strategy.These six strategies are de noted 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  Fig. 2. Compared to the FE analysis, the macroscopic predictions of the APHR method are very good whatever the adopted strategy (the accu racy 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 correla tion 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 accu mulated 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 ele ments.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 corre sponding to the second component of the deformation field ob tained from the APHR model using the REF1 strategy (six modes) and the REF2 strategy (seven modes).In Fig. 6, we can find the dif ferent 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 re duced 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.

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.)

Fig. 1 .
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.)

Fig. 3 .Fig. 4 .
Fig.3.Average acumulated 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 1 .
1 and I dev I I vol so that for a ij = a ji we have: (1)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 u u ve ; vp ; p and r such that:

Table 1
Constitutive model parameters for polycarbonate at 22 °C.

Table 3
Computational times.