Sébastien Mercier
email: sebastien.mercier@univ-lorraine.fr

Katarzyna Kowalczyk-Gajewska

Christophe Czarnota

Effective behavior of composites with combined kinematic and isotropic hardening based on additive tangent Mori-Tanaka scheme

Keywords: elastic-viscoplasticity, homogenization, finite element, metal matrix composite, Mori-Tanaka scheme, kinematic hardening

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

Introduction

Composite materials are widely used in industries. They offer the possibility to combine different materials in order to tailor mechanical properties that a homogeneous material could hardly provide.

Those materials can be manufactured by different processing routes such as powder metallurgy or casting... The development of additive manufacturing could also help providing more dedicated and controled microstructures. In the present work, we focus on two phase elastic-viscoplastic composite.

As a direct application of the contribution, particulate metal matrix composites (which have been widely investigated in the literature) are concerned.

During cyclic loading, heterogeneous materials in general and metal matrix composites in particular may present a significant Bauschinger effect, notably when considering high strength or dual phase steels. Such effect results in the development of kinematic hardening. In the last decades a large literature has been devoted to model such cyclic response based on phenomenological theory, see for instance [START_REF] Chaboche | Mechanics of Solid Materials[END_REF]]. Kinematic hardening observed at the macro level of heterogeneous material can have different sources. Mainly it is a result of the misfit between the mechanical behavior of the constituents. It can also originate in metallic alloys from various dislocation mechanisms observed at the level of a grain (e.g. pile-up dislocations at grain boundaries and around precipitates, [START_REF] Allain | Microstructure based modeling for the mechanical behavior of ferrite-pearlite steels suitable to capture isotropic and kinematic hardening[END_REF]). Capturing the overall response of such materials from the local behavior of the phases has been an issue for many years and is still an open problem. In the present work, we concentrate on the development of a micromechanical scale transition scheme which should include kinematic hardening to the description of the constitutive behavior of the phases.

Estimation of the overall properties of heterogeneous media is a long standing problem in mechanics that, even in the context of linear elasticity has been tackled in numerous ways, see discussion in [Milton(2002), [START_REF] Qu | Fundamentals of Micromechanics of Solids[END_REF], Bruno et al.(2010)]. For discussion on bounds for general composite material, one can refer to [START_REF] Hashin | A variational approach to the theory of the elastic behavior of multiphase materials[END_REF]] (or to [START_REF] Luciano | Non-local constitutive response of a random laminate subjected to configuration-dependent body force[END_REF] ing estimates of the solution to the problem. In the similar spirit, defining the overall behavior of heterogeneous elastic-viscoplastic materials in the real time space is also pursued in [START_REF] Doghri | Mean-field homogenization of elasto-viscoplastic composites based on a general incrementally affine linearization method[END_REF], [START_REF] Agoras | Incremental variational procedure for elasto-viscoplastic composites and application to polymer-and metal-matrix composites reinforced by spheroidal elastic particles[END_REF], [START_REF] Sabiston | Micromechanics based elasto-visco-plastic response of long fibre composites using functionally graded interphases at quasi-static and moderate strain rates[END_REF], [START_REF] Wu | An incremental-secant mean-field homogenization method with second statistical moments for elasto-visco-plastic composite materials[END_REF]]. This way (i.e. solving the problem in the real space) of deriving constitutive behavior for composite materials is followed in the present contribution.

In addition to elastic-viscoplasticity, the incorporation of linear kinematic hardening of the phases into multi-scale approach was recently investigated in the context of variational models [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF], [START_REF] Boudet | An incremental variational formulation for the prediction of the effective work-hardening behavior and field statistics of elasto-(visco)plastic composites[END_REF], [START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF], [START_REF] Lucchetta | A incremental variational procedure for elastoplastic composites with combined isotropic and linear kinematic hardening[END_REF]]. In those formulations the local free energy is described by two potentials, one for the elastic energy and another quadratic one for the stored energy, which is expressed by the (visco)plastic strain tensor and the kinematic hardening moduli. Moreover the dissipation potential has to be modified. Those new terms must be also accounted for in the formulation of the LCC potential which introduces additional difficulties to the minimization procedure.

The aim of this study is to extent the additive tangent Mori-Tanaka scheme proposed for elasticviscoplastic composites to the case where phases exhibit both isotropic and kinematic hardening. The approach combines a tangent linearisation of non-linear viscous response [START_REF] Molinari | A self consistent approach of the large deformation polycrystal viscoplasticity[END_REF]] and an additive interaction law [Molinari(2002), [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]] (or equivalently corresponding sequential scheme [START_REF] Kowalczyk-Gajewska | Sequential linearization method for viscous/elastic heterogeneous materials[END_REF]). Following [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]], the accuracy of the proposed extension is verified with respect to the results of numerical homogenization performed on a representative volume element with randomly distributed spherical particles.

In this paper, we show that from the theoretical point of view accounting for linear or non-linear kinematic hardening within the additive tangent interaction law is straightforward. As in [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF]], only linear kinematic hardening is considered for validation. It will be shown that the resulting micromechanical framework delivers predictions of satisfactory agreement with re-spect to finite element results. The paper is constructed as follows. In section 2 the basic equations of the Mori-Tanaka model based on the additive tangent interaction law and its extension to the case of kinematic hardening are presented. Section 3 includes the methodology and assumptions of the performed finite element calculations. In section 4 which corresponds to the core of the paper, the validation of the proposed micromechanical framework with respect to the results of numerical homogenization is presented. The cases of two-phase composite with hard and soft inclusions subjected to cyclic tension-compression loading are studied. The composite material contains a 25% volume fraction of spherical inclusions.

The present paper is a first step on the way to propose a multi-scheme approach dedicated to heterogeneous elastic-viscoplastic materials with presence of kinematic hardening. As mentioned above, the challenge is quite important so, we first consider a class of metal matrix composites where the matrix phase and the filler are both presenting linear kinematic hardening. In that configuration, it is also shown that the model can be used for elastic-viscoplastic matrix reinforced by elastic inclusions, corresponding to an important class of materials (MMC material). But, for comparison to real materials (see [START_REF] Guo | A cyclic visco-plastic constitutive model for timedependent ratchetting of particle-reinforced metal matrix composites[END_REF]] for T6-treated 6061 Al alloy where the matrix is elastic-viscoplastic and reinforced by SiC particles), non linear kinematic hardening rule must be considered. Application to real material could be investigated in a near future, as a continuation of the present work. It is also worth to note that in most of the recent papers cited above (i.e. [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF], [START_REF] Boudet | An incremental variational formulation for the prediction of the effective work-hardening behavior and field statistics of elasto-(visco)plastic composites[END_REF], [START_REF] Cheng | A homogenization model for porous ductile solids under cyclic loads comprising a matrix with isotropic and linear kinematic hardening[END_REF], [START_REF] Lucchetta | A incremental variational procedure for elastoplastic composites with combined isotropic and linear kinematic hardening[END_REF]]) it is mentioned that the developed approaches are valid for elastic-viscoplastic phases with both isotropic and kinematic hardening. Nevertheless, to the author's best knowledge, no results combining elastic-viscoplasticity and kinematic hardening were deeply investigated in the literature. Thus finite element results included in the paper could also be used for validation of other strategies which exist in the literature.

2 The Mori-Tanaka scheme for composite materials including kinematic hardening

Position of the problem

In the present paper, we consider a two-phase composite material. However, the theory is valid for any composite material with N inclusion phases. The phases have non-linear elastic-viscoplastic behavior, presenting both kinematic and isotropic hardening. Figure 1 The composite material is subjected to a macroscopic strain rate tensor Ė prescribed at the remote boundary of the RVE. The total strain rate in each phase is decomposed into elastic and viscous parts:

ε = εe + εv . (1)
The elastic strain rate εe is linked to the Cauchy stress rate tensor σ by the incremental elastic law:

εe = M e • σ , (2)
where M e is the fourth order tensor of elastic compliance and the operator • is the double contracted product. The present contribution is restricted to a small deformation theory. The viscoplastic contribution εv is assumed volume preserving. At the level of each phase, a constitutive law accounting for viscoplasticity and kinematic hardening is considered. In this contribution, Perzyna-type viscoplasticity model with an over-stress function Φ is adopted, [Perzyna(1986)] :

εv = Φ ∂f ∂σ = Φ 3(s -x) 2σ eq , (3)
where s represents the deviator of the Cauchy stress tensor, x is the backstress tensor and

f = σ eq -σ Y -R(ε eq ) , Φ =                εo f σ Y +R(ε eq ) 1 m if f > 0, 0 otherwise . (4)
Since kinematic hardening is accounted for in the present work, the equivalent Huber-von Mises stress is defined as σ eq = 3/2(sx) • (sx). The function R(ε eq ) describes the hardening behavior of the phase. Variable ε eq = Φdt is the accumulated plastic strain where Φ = εeq = 2/3 εv • εv is the equivalent viscoplastic strain rate. For this model, the viscous flow is initiated when the function f is positive. Parameters m and εo are representative of the rate dependency of the flow law, σ Y is the initial flow stress under static conditions. The function R(ε eq ) is taken as non-linear :

R(ε eq ) = k(ε eq ) n , (5)
where k and n are material parameters. From Eqs (3) to ( 5), the Perzyna model can be rewritten in the following form (valid after initiation of the viscous flow):

σ eq = (σ Y + k(ε eq ) n ) 1 + εeq εo m . (6)
As proposed by [Prager(1949)], a linear kinematic hardening law is adopted in the present paper :

ẋ = 2 3 h εv (7)
in which h is the linear kinematic hardening modulus. Nevertheless, since the proposed Mori-Tanaka scheme is also valid for non linear kinematic hardening, the kinematic law of [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF]] (see also [START_REF] Frederick | A mathematical representation of the multiaxial bauschinger effect[END_REF]]) or [Chaboche(1989)] can be assumed as well.

We focus on the modeling of composites reinforced with spherical inclusions and propose to adopt the Mori-Tanaka scheme based on [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]] (see also [START_REF] Kowalczyk-Gajewska | Sequential linearization method for viscous/elastic heterogeneous materials[END_REF]]) to account for kinematic hardening. The model is based on the tangent linearization of the viscous response of the phases, see [START_REF] Molinari | A self consistent approach of the large deformation polycrystal viscoplasticity[END_REF]]. As a consequence, it has to be checked if the tangent linearization of the phase behavior is still preserved in its mathematical form when kinematic hardening is added. This will be shown in the next section.

Additive tangent interaction law with kinematic hardening

When plasticity is activated, the viscous strain rate tensor can be written in the following form, based on Eqs. ( 3) and (4) :

εv = 3 εo 2σ eq σ eq -σY σY 1 m (s -x) , where σY = σ Y + R(ε eq ) (8)
The fourth order tensor for viscous tangent compliance is defined as :

M v(tan) = ∂ εv ∂σ = ∂ εv (σ, x, ε eq ) ∂σ (9)
From Eq. ( 8), one obtains :

M v(tan) = 3 εo 2σ Y σ eq -σY σY 1 m -1 3 2 1 m - σ eq -σY σ eq s -x σ eq ⊗ s -x σ eq + σ eq -σY σ eq K ( 10
)
where K is the fourth order identity tensor for deviatoric second order tensors s, i.e. K • s = s and K • σ = s, ⊗ is the tensorial product operator. Note that the definition ( 10) is identical to the expression of the inverse tangent stiffness for a phase without kinematic hardening when (sx) is replaced by s. In addition, the kinematic backstress is accounted for in this relationship through the expression of the equivalent Huber-von Mises stress.

From the definition of the viscous tangent compliance tensor, the term M v(tan) .(sx) can be evaluated :

M v(tan) .(s -x) = 1 m σ eq σ eq -σY εv (11)
Therefore, the linearized form of viscoplastic relation ( 8) is obtained:

εv = M v(tan) .σ + εv res ( 12
)
with

εv res = m σ eq -σY σ eq -1 M v(tan) .(s -x) -M v(tan) .x (13)
is a back extrapolated strain. Following the standard procedure of mean-field homogenization, for further specification of a scale-transition scheme which relays on the Eshelby result, it is assumed that M v(tan) and εv res are calculated for averaged stress and backstress in the phases.

We have shown that the introduction of the backstress in the constitutive behavior still enables to define a tangent linearization which is consistent with the approach proposed by [START_REF] Molinari | A self consistent approach of the large deformation polycrystal viscoplasticity[END_REF]].

As a consequence, the Mori-Tanaka scheme based on the additive tangent interaction law and proposed firstly in [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]] can be adopted.

Adopting the same reasoning as for the elastic-viscoplastic materials without kinematic hardening, with the additive tangent interaction law, the average strain rate εi in the inclusion can be linked to the average strain rate εm in the matrix, see Fig. 1, as follows :

εi -εm = -M v * • (s i -s m ) -M e * • ( σi -σm ) , (14)
σ m (resp. σ i ) represents the average Cauchy stress tensor in the matrix (resp. in the inclusion). s m (resp. s i ) denotes its deviatoric part. M v * and M e * are the inverse Hill tensor obtained with use of the viscous tangent and elastic stiffness. Note that when the first term of the right hand side of Eq. ( 14) related to the viscoplastic tangent contribution is disregarded, the Kroner model is retrieved.

Therefore, in Section 4, the difference in predictions between our new approach and the Kroner model will be discussed in the frame of kinematic hardening contribution.

From the interaction law ( 14) and from the incremental elastic law (2) for inclusion and matrix phases, the average strain rate in the inclusion is given by:

εi = (M e -1 i + M e -1 * ) -1 • M e -1 i • εv i + M e -1 m • ( εm -εv m ) + M e -1 * • { εm -M v * • (s i -s m )} (15)
From the consistency equation < ε >= Ė, the strain rate in the matrix εm is found and henceforth, the strain rate in the inclusion εi is obtained. The macroscopic stress tensor Σ and the corresponding deviator S are obtained by volume averaging over the whole RVE:

Σ =< σ > , S =< s > . ( 16
)
More details can be found in [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]] and in an equivalent formulation of the Mori-Tanaka scheme from the sequential approach proposed in [START_REF] Kowalczyk-Gajewska | Sequential linearization method for viscous/elastic heterogeneous materials[END_REF]].

The time integration procedure of the proposed Mori-Tanaka scheme follows a forward Euler method.

The main steps are summarized in Appendix A for a two phase composite subjected to cyclic tension/compression. Note that as mentioned before, the present approach can be extended to a composite with N inclusion phases with no additional effort, since Eqs ( 14) to (16) remain valid. In addition, the present development could also be used for the development of a self consistent scheme for polycrystalline materials owing to the fact that the additive tangent interaction law ( 14) is preserved when accounting for kinematic hardening.

Numerical model

In the previous section, it has been shown that the Mori-Tanaka scheme based on the additive tangent interaction law can be extended to the case of two phase materials having elastic-viscoplastic behavior with combined kinematic and isotropic hardening. The extension is straightforward and does not involve large modification in the numerical implementation of the model. Nevertheless, the proposition has to be validated. This step will be carried out having recourse to finite element simulations with Abaqus Software. The case of a composite material with c = 0.25 volume fraction of inclusions is considered.

The 3D representative volume element with 30 inclusions, as proposed in [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]] is used. The statistical homogeneity and isotropy of the adopted RVE has been verified using two point probability functions, see [START_REF] Mortazavi | Modeling of two-phase random composite materials by finite element, Mori-Tanaka and strong contrast methods[END_REF]] or [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]]. Note also that such a RVE with 30 inclusions is often adopted in the literature to validate predictions of multi-scale models, see for instance [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF]] and [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF]]. The RVE is meshed with 10-node full integration quadratic tetrahedral elements (C3D10 in ABAQUS). The mean number of elements per inclusion is 1, 250 and the whole cell contains 150, 000 elements, see Fig. 2a. It has been shown in [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]] that such mesh density ensures mesh independent results for the average values in the phases and in each inclusion. Periodic boundary conditions are also imposed at the external surfaces of the RVE for the present work, following the numerical technique proposed by [Al Kassem(2010)] for non conformed mesh. Note that periodic boundary conditions have also been considered elsewhere, see [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF]]. Geometric nonlinearities are disregarded since a small deformation formalism has been adopted in the theory.

The elastic-viscoplastic Perzyna model with kinematic hardening is not implemented in Abaqus Standard. As a consequence, a UMAT subroutine has been developed, following the strategy proposed by [START_REF] Lubarda | On the numerical algorithm for isotropickinematic hardening with the Armstrong-Frederick evolution of the back stress[END_REF]] or [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]]. The main steps of the integration procedure are recalled in Appendix B. In the present paper, only linear kinematic hardening is con- sidered. Validation for non linear kinematic hardening could be tackled in a future extension of the present work.

Results

As mentioned previously, the MT model has been developed to describe the overall behavior of the composite response. Our goal is to show that for a quite large volume fraction of inclusions (here c = 0.25), the predictions remain accurate at the macroscopic scale but also for the mean local response of the phases. The composite material is subjected to tension/ compression tests. The longitudinal strain rate is 10 -3 s -1 . The strain amplitude is fixed to E max 11 = 0.01. A larger strain amplitude E max 11 = 0.04 was also tested in [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]] where it was observed that the predictions of the model was more accurate for large strain because plasticity was able to develop in the entire RVE. With a smaller strain amplitude, a confined plasticity may develop depending on the respective material properties.

So adopting this small strain amplitude is more challenging for the model than considering larger strain.

In the following, we present results for stress, backstress and strain components in the direction of the loading 1. They correspond to average value at the level of each individual phase or of the composite.

In the following of the paper, when a mention is made to stress, backstress or strain, this will always be referring to the above components, if not otherwise stated.

Two phase composite with hard inclusions

The reference material parameters are given in Table (1). The inclusion phase is hard. The same moduli for the linear kinematic law (7) are adopted for the matrix and inclusion phases, h = 10GP a.

Table 1 Case of hard inclusions. Material parameters for the inclusion and matrix phases. Viscoplasticity is described by the Perzyna-type law (6). Most of the parameters are taken from [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF]].

Phase To evaluate the quality of the predictions also at the inclusion level, the mean stress versus strain response in the inclusion is additionally displayed in Fig. 3 b). From the FE comparisons, one observes that for the first loading/unloading, the predictions of our additive tangent MT model is quite consistent with the FE mean value. Indeed, for the first loading, the predictions of the Finite Element model shows that limited plasticity has developed. With our model, the response remains elastic while for the Kroner model, the mean value of the strain in the inclusion domain reaches almost 0.004. Clearly, the plastic strain accumulated during the first loading is too large with the Kroner approach when compared to the FE results. Note also that the local stress in the inclusion phase is well captured by our model from the first to the last cycles, while the Kroner predictions show much more discrepancies. As the cyclic loading is continued, the difference between the two MT schemes and FE results is reduced.

E [GPa] ν σ Y [MPa] k [GPa] n ε0 [s -1 ] m h [GPa]
To have a better understanding of the evolution of the mechanical fields during cycles, the time evolution of the longitudinal stress and backstress components in the two phases is presented in Fig. 4. From finite element calculations, the mean backstress is calculated by volume averaging over the elements belonging to the domain of interest (here inclusion or matrix). For the matrix phase, the stress evaluation shown in Fig. 4 a) is correct with the Mori-Tanaka approach based on the additive tangent interaction law and the Kroner based approach, when compared to the FE calculations. We also observe that the backstress component x m 11 in the matrix phase presented in Fig. 4 b) is accuratly predicted by our approach from the first cycle, while for the Kroner based approach, the predictions become accurate only after the first loading/ unloading. Since a linear kinematic hardening law (5) is accounted for in the paper, the backstress evolution is an indicator of the discrepancy (or accuracy)

regarding the level of the plastic strain rate. From 1.

development calculated by the additive tangent MT scheme is too small at the early stage of the cyclic loading (almost elastic response during the first cycle), while in the FE model, plastic activity starts from the first cycle. This delay in the initiation of the plastic deformation observed with our model has an effect on the following cycles. Nevertheless, as the number of cycles is increased, the observed difference between our multi-scale approach and FE results shrinks, meaning that the predictions in terms of plasticity development, and as a consequence, of the mean strain and stress, are evolving to become accurate. With the Kroner approach, the mean stress component σ i 11 and backtress component

x i 11 in the inclusion domain are too large during the early deformation stage when compared to the FE results. This is due to a larger accumulated plastic strain. This effect is already depicted in Fig. 3 b)

where the strain amplitude in the inclusion domain during the first loading is by far too large. After two cycles, the prediction of the Kroner based Mori-Tanaka scheme seems to be quite consistent. Still, our approach provides more salient predictions.

It is also observed in FE calculations that the magnitude of the matrix backstress decreases as the number of cycles increases, while in the inclusion it first increases and then decreases after the second cycle. Remember that the fluctuation of the backstress is directly proportional to the viscoplastic strain rate. Indeed, after few cycles, due to isotropic hardening, the intensity of viscoplastic flow and its time duration within a cycle decrease since the RVE is subjected to constant cyclic strain amplitude. Because of the isotropic hardening law which is not saturating with strain (see Eq. ( 5)), it is expected that after a large number of cycles, the response of the two phases becomes elastic, see Figs

b) and d).

As in [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]], FE results obtained in this paper enable to capture heterogeneity developed in the inclusion domain. For a given field variable (strain, stress, backstress tensor components), a first level of heterogeneity is obtained by considering the mean value of the variable in each of the 30 inclusions. A second level of heterogeneity portrays probability density functions at the finite element scale, giving a clear indication on how the variable fluctuates in the inclusion domain, see Appendix D. Results derived by the proposed MT scheme are also provided together with the mean value evaluated from the numerical simulation. The heterogenity in the stress field appears quite limited for all cycles where a good match is observed between the MT approach and the FE results, see Fig. Indeed, in the Finite Element calculations, it is observed (see also [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]]) that during the first cycle, since the overall prescribed strain amplitude is quite small, plasticity will develop in a heterogeneous manner within the inclusion domain. This type of response is really difficult to capture by mean field approaches, explaining why our MT approach suffers in restituting the situation at the 1.

early stage of the cyclic loading. Nevertheless, after some cycles, heterogeneity is strongly reduced and the MT approach reproduces quite satisfactorily FE strain and backstress, as already depicted in The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1. exhibits a normal like distribution shape around the mean value, see Fig. 6 b). In contrast, Fig. 6 c)

shows for the backstress that the fluctuation between inclusions is particularly marked and the corresponding probability density function at the finite element level is of exponential type, see Fig. 6d). Note that further during the loading process, when plasticity has developed in all inclusions, the probability density function in backstress will evolve from an exponential type to a bell shaped one. A careful inspection of Fig. 6c) enlights the fact that 3 inclusions remain almost elastic, 16 inclusions present a value of the average backstress which is quite lower than the average value (see solid line), 3 inclusions are presenting value close the average one while only 8 develop more plasticity than the average. For inclusions presenting limited backstress (3 inclusions here behaving almost elastically), the average strain in such inclusion is close to the one predicted by the MT approach since for our Mori-Tanaka model, plasticity is not triggered during the first loading, see Fig. 6e).

Various measures are next introduced to quantify more precisely the evolution of the heterogeneity in FE calculations. Let us introduce r σ , r ε , r x as measures of absolute deviations around the mean value at the finite element level and r I σ , r I ε , r I x the corresponding measures at the level of the inclusion phase, see Eqs (D.4-D.5) of Appendix D for their definition. In addition, the absolute relative difference in inclusion mean longitudinal strain between our MT approach and the FE calculations is defined as:

d ε 11 = |ε MT 11 -ε FEM 11 | 1 2 |ε MT 11 + ε FEM 11 | (17)
where ε MT 11 (resp. ε FEM 11 ) stands for the average longitudinal strain in the inclusion phase predicted by the additive tangent MT approach (resp. by the FE calculation). Using a definition analogous to (17), d σ 11 and d x 11 are also considered in the following.

Fig. 7 shows absolute deviations in terms of longitudinal stress (a), backstress (b) and strain (c) at the two levels of observation taken at each point of the loading reversal. For comparison, the absolute difference (17) between the additive tangent MT approach and the FE results is also displayed. Fig. 7 notably shows that relative absolute deviations at the inclusion level, r I σ , r I x and r I ε are lower than the corresponding r σ , r x and r ε at the level of finite elements. Narrower fluctuations between inclusions are due to the fact that the heterogeneity in stress, backstress and strain, observed at the finite element level (see Fig. 6a) is partially reduced when averaging at the inclusion level . However, the measures at the two scales of heterogeneity reveal very similar trends. When concentrating on the stress, Fig. 7a) depicts the reduction in the fluctuations at the inclusion and finite element scales as time is increased with asymptotic values of 2.6% for r I σ and 5.8% for r σ . Correspondingly, the absolute difference d σ 11

between our mean-field approach and the average FE value in the inclusion phase is lower than 2.5% throughout the entire cyclic loading. This illustrates again the great capability of the proposed modeling to describe the stress in the inclusion phase, for the case of hard inclusions adressed in this section, see also Fig. 3b). Fig. 7b) and c) related to the backstress and strain respectively, reveals that when the heterogeneity is particularly marked in FE (for t = 10, 30s), the absolute difference between FE results and MT approach is larger too. Nevertheless, from Fig. 7a) to c) no direct correlation between the heterogeneity and the absolute difference is clearly highlighted. As an illustration, from Fig. 7a), one observes that d σ 11 = 0.2%, r I σ = 2.7% and r σ = 6.4% at t = 50s, while at t = 130s d σ 11 is increased up to 2.3% with a slight decrease of both r I σ = 2.6% and r σ = 5.7%. This may suggest that the mismatch is not only due to the heterogeneity developed at various scales in the FE model. The delay in the development of plastic deformation in the inclusion phase observed in the MT approach, see Fig. 5-c), may be partly at the origin of the discrepancy between FE results and our modeling.

For sake of brevity, the predictions without kinematic hardening are not displayed here since they can be found in Fig. 11 of [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]]. It has to be mentioned that for the present set of material parameters (most of the parameters are taken from [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF]]), better predictions are obtained with the adopted value of linear kinematic hardening parameter, especially for the strain development in the inclusion phase.

A good quality of the predictions after some cycles, as demonstrated in Figures 3 and4, could be partially attributed to the plastic accommodation. Of course, hardening plays a positive role in rendering the response more homogeneous, and therefore, predictions by mean field approaches eased.

To understand in a deeper manner the influence of isotropic hardening, lower values of the parameter k involved in the description of the isotropic hardening of both phases are considered, namely k = 2.5GP a and k = 0GP a. For the latter case, no isotropic hardening is present and the flow stress is simply rate dependent of the form: 1 except for the isotropic hardening coefficient involved in Eq. ( 5)

σ eq = σ Y 1 + εeq εo m . (18)
whose value is k = 2.5GP a.

For k = 2.5GP a, it is clearly seen in Fig. 8 that the overall stress-strain response under cyclic loading is still well predicted. The transition from an elastic to a viscoplastic response is also well reproduced.

We have to mention that the difference between FE and MT is decreasing as the number of cycles is increased. At the level of the inclusion, the strain amplitude during a cycle becomes well captured after 3 cycles. This is quite equivalent to the previous case where the isotropic hardening coefficient was set to k = 5GP a. Nevertheless, one sees that the difference between FE and MT is enlarged as the isotropic hardening is reduced. For the time evolution of the mean stress and backstress in the matrix, with k = 2.5GP a, no significant mismatch is observed in Fig. 9 The lack of isotropic hardening induces more heterogeneity in the strain field leading to part of the domain which remains elastic and part of the domain which has already plastified. As a consequence, it is more difficult for any mean field approach to capture this stage. The Kroner based model with kinematic hardening is not able to provide consistent results, as seen in Fig. 10. For the inclusion response, the additive tangent MT scheme predicts that the inclusion domain remains elastic while in the FE calculations, it is observed that some (very limited) plasticity may develop in inclusions.

For a Kroner based model, the strain amplitude in the inclusion domain is too large during a cycle.

Clearly, the inclusion phase develops large plasticity which is not observed in the FE calculations.

This situation explains the poor prediction capabilities of the overall response of Fig. 10 inherited from the Kroner based approach.

Next, we concentrate on the time evolution of the stress and backstress within the inclusion domain, since the time evolution of the stress and backstress in the matrix phase is accurate with the proposed additive tangent MT scheme (results not presented here). Fig. 11 a) shows that the stress level is almost accurately predicted with our model while with the Kroner approach, a too stiff response is obtained. Concerning the backstress evolution provided in Fig. 11 b), a non-zero mean value of the finite element backstress in the particles exits (less than 1 MPa in magnitude), while the MT prediction is null, meaning that no plasticity is triggered. Due to a large strain value predicted in the inclusion domain by the Kroner based approach, the backstress is too large (maximum value of the order of 30 MPa).

When the contribution of the kinematic hardening is also disregarded (no hardening), the two phases have an elastic ideally-viscoplastic response. It is observed in Fig. 12, that since the inclusion is the hard phase with an initial yield stress five times larger than the matrix phase, no plasticity will be experienced by the inclusion phase. So the response is purely elastic. For the matrix and overall responses, the predictions are still correct, with some mismatch during the elastic-viscoplastic transition. It is also observed that even for a cyclic loading of low strain amplitude E max 11 = 0.01, the plateau regime corresponding to the steady state for viscoplasticity where elastic contribution vanishes, is reached. In that case, the good correspondence between the additive tangent MT and the FE models confirms that the interaction law is accurate for viscoplastic flow with hard inclusions, see [START_REF] Molinari | A self consistent approach of the large deformation polycrystal viscoplasticity[END_REF]], [START_REF] Molinari | Tuning a self-consistent viscoplastic by finite elements results -I modeling[END_REF]] or [START_REF] Molinari | Validation of the tangent formulation for the solution of the non-linear eshelby inclusion problem[END_REF]].

The role of kinematic hardening is further investigated by prescribing larger values of the kinematic The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1 except for the isotropic hardening parameter k which is set to k = 0GP a.

stress-strain response based on the additive tangent MT scheme are compared to FE results for the three values of h and presented in Fig. 13. Only the response during the first cycle is presented here.

Since in the previous calculations of this section, it was observed that the less accurate predictions were obtained during the first cycle, this presentation is fair to judge whether the model is able to precisely predict the overall response for large linear kinematic hardening parameter. It is observed that as h is increased, our model is still able to reproduce the stress-strain response with a good accuracy.

The time evolution of the stress in the inclusion phase is next presented in Fig. 14 is harder than the matrix, with material parameters given in Table 1. The kinematic hardening modulus is varied from h = 10GP a to h = 10 3 GP a. The volume fraction of inclusions is c = 0.25.

domain is also in a good agreement with the FE results. In addition, it is also observed that as the kinematic hardening modulus is becoming large, the role of isotropic hardening is reduced. For the most extreme case with h = 1000 GPa, we can observe that the time evolution of the longitudinal stress is almost linear. In Fig. 14 b), the local backstress in the inclusion phase is also depicted. For all three cases, we observe that the predictions are becoming more accurate as the number of cycles increases. For the fifth cycle (which is the last one considered in this calculation), the agreement is almost perfect. Note that through the backstress evolution, because of linear kinematic hardening, this figure proves that the plastic activity for the inclusion phase, captured in average by the proposed Mori-Tanaka scheme, is also in a good agreement with the FE results. It must also be mentioned that the Kroner based model and the additive tangent model are providing similar results for large h in terms of overall response. In Appendix C, the Kroner based Mori-Tanaka scheme with kinematic hardening and without isotropic hardening is investigated and the response under uniaxial tension is derived in an approximated way. It is observed that as h is increasing, the slope of the macroscopic response of the composite is tending to an elastic one, meaning that when h is large, plasticity vanishes and a linear response is recovered. This finding explains why as h is large, the hysteresis loop during cycle is reduced and tends to vanish, see Fig. 13.

The last case investigated in this section dedicated to hard inclusion phase is the limit case of elastic inclusions. The matrix phase has an elastic-viscoplastic response with combined isotropic and kinematic hardening. Material parameters are those presented in Table 1, except for the inclusion which remains elastic. This case is interesting since it corresponds to the configuration of a metal matrix composite with elastic inclusions (for instance SiC particles). Nevertheless, for a real material, the kinematic hardening rule is non linear. In [START_REF] Guo | A cyclic visco-plastic constitutive model for timedependent ratchetting of particle-reinforced metal matrix composites[END_REF]], a model where the total backstress is divided in ten parts (the evolution law being non linear), is adopted to reproduce accurately the behavior of the matrix phase (T6 -treated 6061 Al alloy). Such extension to non linear kinematic rule will be considered in a near future. The overall response and the stress-strain behavior in the inclusion phase are provided in Fig. 15 a). The predictions are accurate. The matrix response is presented in Fig. 15 b). It is seen that the time evolution of the predicted stress during the cycles fits the FE results.

In addition, the mean backstress in the matrix phase evaluated by the proposed mean field and by the MT h=10 GPa FEM, h=10 GPa MT h=100 GPa FEM, h=100 GPa MT h=1000 GPa FEM, h=1000 GPa Fig. 14. Time evolution under cyclic loading of the longitudinal a) stress and b) backstress in the inclusion phase. The prediction of the proposed MT scheme is compared to the FE results. The inclusion is harder than the matrix, with material parameters given in Table 1. Three values of the kinematic hardening modulus h from h = 10GP a to h = 10 3 GP a are considered. The volume fraction of inclusions is c = 0.25.

plastic strain rate and also the strain are accurately predicted in the matrix phase at least in average.

Note that for subsequent cycles, due to the presence of isotropic hardening, the time interval where the matrix remains elastic increases. For a large number of cycles, the matrix response will approach a purely elastic one. That is also why a better agreement between FE and MT predictions is observed with increasing the number of cycles. sions is c = 0.25. The material parameters are provided in Table 1 except for the inclusion phase which is purely elastic.

Two phase composite with soft inclusions

In this section, the case of soft inclusions is investigated. The material parameters are provided in Table 2 and are similar to those of Table 1 with switching matrix/inclusion properties. So the contrast in the Young moduli and for the initial yield stress is 1/5. As it is seen in Fig. 16 for the composite response under the applied cyclic loading, significant discrepancy exists mostly in terms of stress amplitude. Interestingly, the strain in the inclusion domain is well captured but the stress level is not so satisfactory. One can mention that in the case of hard inclusions, it is often the stress level which is Case of soft inclusions. Material parameters for the inclusion and matrix phases. Viscoplasticity is described by the Perzyna-type law (6).

Phase well captured while the strain amplitude in the inclusion phase is not so precisely predicted, as seen in the previous section 4.1, see also [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]]. In Fig. 17 a) for the matrix phase, the estimates of stress and backstress are really satisfactory. As time increases, the stress level increases mostly due to the hardening. In addition, the magnitude of the backstress is slowly reduced, showing that less plasticity is accumulated as the number of cycle is getting larger. For the corresponding quantities in the inclusion phase on Fig. 17 b), the stress level is underestimated by our multi scale approach. Nevertheless, the predictions become more and more accurate during subsequent cycles, again due to isotropic hardening. The backstress levels are quite similar for the additive tangent Mori-Tanaka scheme and FE calculations. The reduction of the backstress magnitude implies that the magnitude of the viscoplastic strain is reduced. While initially, the component of viscoplastic strain is larger in the inclusion phase, becomes more important in the matrix phase after three cycles. In addition, the duration where the plastic flow is present in the inclusion is reduced (direct effect of the lower Young modulus in the inclusion phase). It could be anticipated that the inclusion will remain mostly elastic in few subsequent cycles. Note that while the initial yield stress ratio is 1/5 and due to considering the same isotropic hardening parameters for both phases (same n and k as seen in Table 2), the ratio of flow stress between phases is reduced with cycles, reaching a value about 1/3 at the end of the last tensile ramp.

E [GPa] ν σ Y [MPa] k [GPa] n ε0 [s -1 ] m h [GPa]
We now analyze the heterogeneity developed in the FE model for the case of soft inclusions. 2. The inclusion phase is soft.

value evaluated from the FE calculation and from the proposed additive tangent Mori-Tanaka scheme is also added. It clearly appears from Fig. 18a), that the fluctuation in stress from one inclusion to another is very limited during the whole deformation process. Fig. 18a) also confirms that the MT approach has difficulties in terms of stress evaluation to capture precisely the elastic-plastic transition for the soft inclusion configuration, at least with the adopted material parameters. This is particularly marked at the beginning of the cyclic loading, see Let us consider now the variables introduced previously: d κ , with κ = {σ i 11 , x i 11 , ε i 11 }, related to the accuracy of the MT approach against FE results (Eq. 17 As already mentioned, the measure at the finite element level presents a larger deviation when compared to the deviation evaluated based on the 30 inclusions. Concentrating now more specifically on the longitudinal stress, it appears from Fig. 20a) that r I σ and r σ are nearly constant from the first cycles with mean values, calculated over the entire process, around 3% for r I σ and 5% for r σ . Heterogeneities in stress at the inclusion or at the finite element scales thus are very limited even at the beginning of the process. In contrast, the absolute difference between the MT approach and the FE calculation is decreasing with d σ 11 ≃ 30% at t = 10s and d σ 11 ≃ 6% at the end of the process. The situation for the backstress depicted in Fig. 20b) is somewhat different. Relatively low values of about 3% for r I

x and 6% for r x are determined over the entire process. At the late stage, the relative absolute deviations around the mean is seen to increase. This can be explained as follows. As time increases, it is observed in Fig. 17 b) that the average value of the backstress is reduced. So the increase in r I x and r x is caused by the lower value of the mean backstress. It has been checked that the absolute deviations which can be defined as |x i 11 |r I x and |x i 11 |r x are reaching asymptotic values below 0.5M P a as time is large. For the longitudinal strain, a correspondence may be found in Fig. 20c) between evolutions of the absolute difference d ε 11 and absolute deviations r I ε and r ε . The corresponding values are decreasing when analyzing these indicators at subsequent reversal points up to 70s, then later fluctuating around 2.7%

for r I ε , 5% for r ε and 0.3% for d ε 11 . The results presented in Fig. 20 for the case of soft inclusions (as for the hard inclusion case), suggest that the mismatch between MT and FE results cannot be solely related to the heterogeneity of a given field observed in FE calculations.

Since hardening of the phases may modify the phase contrast (at least in terms of flow stress) during a)

time (s) deviations and difference in stress 2, are those of Fig. 16.

loading, then it has been decided to run similar calculations with k = 0GP a so that no isotropic hardening is included in both phases. As observed for the hard inclusion case, with no isotropic hardening and because of linear kinematic hardening, as expected, the predictions are stabilized after one cycle. Fig. 21 a) presents the macroscopic response of the composite. It is seen that on average, the stress level is well captured at large strain. The difference is observed mostly during the elastic-plastic transition. In our multi-scale approach, the transition is more abrupt, while in the FE calculations, a smooth transition takes place. This is mainly due to the heterogeneity within the phases. Fig. 21 b)

displays the stress-strain response in the inclusion domain. While the strain amplitude is well predicted, the stress level is underestimated by our approach. In the present Mori-Tanaka scheme, the elastic viscoplastic response of the inclusion phase presents three different stages (almost linear by parts). The first linear stage is a signature of a purely elastic response. In a subsequent stage, an affine response is observed when the inclusion phase develops plasticity while the matrix remains elastic.

The third stage corresponds to a situation where both phases are plastified. Fig. 22 displays the time evolution of the stress and backstress in both phases, separately. For both phases, the longitudinal backstress is in general well predicted, meaning that the duration of the plastic flow and intensity are accurately captured. When analyzing carefully the time evolution of the backstress, we can see that the duration of the plateau present for both phases is different. During the plateau stage, the value of the backstress is frozen, so the associated viscoplastic strain rate vanishes. For the longitudinal stress level, Fig. 22 a) shows that the stress in the matrix phase is slightly underestimated while in Fig. 22 b), the stress in the inclusion domain is significantly overestimated by the proposed additive tangent Mori-Tanaka scheme. A supplementary calculation has been performed where the hardening (both kinematic and isotropic) is disregarded: k = h = 0GP a (results not presented here). The goal is to evaluate whether the trend observed in Figs 21 and 22 can be related to the presence of kinematic hardening or not. By analyzing the predictions, we can conclude that for the set of material parameters considered in the present work, the relative mismatch between the theoretical and numerical predictions for these two cases (with or without isotropic hardening) is transferred without any significant changes in terms of accuracy, when kinematic hardening is introduced.

Conclusion

The main goal of the present work is to explore the possibility to extent the additive tangent interaction law proposed initially for elastic-viscoplastic materials to the context of kinematic hardening.

This extension is shown to be straightforward and does not need important changes in the theoretical framework when compared to the previous model dedicated to materials presenting no kinematic hardening, see [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]] or [START_REF] Kowalczyk-Gajewska | Sequential linearization method for viscous/elastic heterogeneous materials[END_REF]]. As a consequence, the numerical implementation of the new model is direct. At this stage, linear kinematic hardening is only considered, as in [START_REF] Lahellec | Effective response and field statistics in elastoplastic and elasto-viscoplastic composites under radial and non-radial loadings[END_REF]]. Extension to non linear kinematic hardening is from a theoretical point of view, possible and could be integrated in a near future.

For two phase composite, a Mori-Tanaka scheme is next formulated and the validation of the proposed approach is performed by carrying out extensive comparison with FE models. For that purpose, the representative volume element developed in [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]] and containing 30 inclusions is considered. In this contribution, periodic boundary conditions are prescribed at the remote boundary of the RVE.

As a reference, a composite with 25% volume fraction of inclusions is subjected to overall uniaxial cyclic tension compression loading, prescribing a small strain amplitude at the remote boundary. For elastic-viscoplastic materials, the most difficult part to be modeled precisely is the transition from elastic to viscoplastic regime. Therefore, a small strain amplitude of E max 11 = 0.01 is adopted here.

The large volume content of inclusions assumed in the present work is also challenging.

The cases of hard and soft phases for the inclusion are considered separately. For the hard inclusion case studied in the present work, the overall response of the composite is accurately described, especially when isotropic hardening is present. When isotropic hardening is disregarded, the agreement between the model and the FE calculations is still consistent. However, the difference between the proposed model and the FE results is larger. Nevertheless, all the key ingredients (strain amplitude, stress level, duration of the plastic loading during cycle, level of the backstress in both phases) displayed by the FE calculations have been captured by the present framework. For soft inclusions, the prediction of the additive tangent Mori-Tanaka scheme is less accurate when isotropic hardening is disregarded. One should recognize that this is a general feature for most of the multi-scale approaches.

Interestingly, some discrepancies were already existing for the same local constitutive laws without kinematic hardening, see [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]]. Therefore, from the present parametric study, it is recognized that the quality of predictions of the additive tangent MT model is not affected (or not induced) by the introduction of the kinematic hardening component in the constitutive behavior of the local phase.

The finite element results provide additional information concerning the heterogeneity which develops in the RVE during loading. Such heterogeneity cannot be captured by the present MT approach.

Nevertheless, it is seen that for the hard inclusion case when isotropic hardening is introduced, the scatter in stress and strain fields is limited and the mean value is really representative of the situation within the material. For the soft inclusion case, the prediction of the model in terms of stress level within the inclusion phase is not so accurate. By defining absolute deviation indicators for three mechanical fields (stress, strain and backstress), it is observed that the mismatch between the proposed Mori-Tanaka and the mean value obtained from FE simulations is not always due to the heterogeneity which prevails inside the RVE.

Even if the predictions of the proposed model are shown to be in general consistent with respect to FE model, some differences still exist which may originate at least from two sources. First, the Mori-Tanaka scheme is known to be salient for composite material with inclusions of moderate volume content, while here the considered volume fraction is quite large. This choice is done on purpose since our goal is to perform a validation of the proposed model, in a demanding configuration. Second, the difference may also come from the interaction law which is proposing an estimate of the Eshelby solution. Therefore, it is of interest to investigate also the accuracy of the interaction law, in the spirit of [START_REF] Mercier | Validation of an interaction law for the eshelby inclusion problem in elasto-viscoplasticity[END_REF]] for elastic-viscoplastic materials, when kinematic hardening is present. This work is under progress.

The Cauchy stress tensor at time t + dt can be split in a deviatoric part and a spherical part:

σ t+dt = s t+dt + 1 3 tr(σ t+dt )I (B.1)
where I is the second order identity tensor. The increment of the total strain can also be decomposed into a deviatoric part ∆e and a spherical part tr(∆ε). The update of the spherical part of the Cauchy stress tensor is straightforward (no volumetric plasticity) so tr(σ t+dt ) = Ktr(∆ε) while the deviatoric part of the Cauchy stress tensor at time t + dt is expressed as:

s t+dt = s t + 2G(∆e -∆ε v ) (B.2)
From Eq. ( 7), the increment of the backstress is proportional to the increment of the plastic strain tensor ∆ε v , so that:

x t+dt = x t + 2h∆ε v 3 (B.3)
∆e is known while the increment of the viscoplastic strain tensor ∆ε v is evaluated via the radial return method. From the definition of the elastic predictors for the deviatoric Cauchy stress tensor s pr t+dt and the backstress tensor x pr t+dt : s pr t+dt = s t + 2G∆e and x pr t+dt = x t , the update of the deviator of the Cauchy stress tensor and of the backstress tensor are only related to the increment of the viscoplastic strain tensor ∆ε v which is the only unknown quantity:

s t+dt = s pr t+dt -2G∆ε v x t+dt = x pr t+dt + 2h∆ε v 3 (B.4)
Besides, we define the known tensor B pr as:

B pr = s pr t+dt -x pr t+dt = s t -x t + 2G∆e (B.5)
From Eqs (3) and ( 4), the increment of the viscoplastic strain tensor is:

∆ε v = 3 2 s t+dt -x t+dt σ eq t+dt Φdt with Φ = εo σ eq t+dt -σ Y -R(ε eq t+dt ) σ Y + R(ε eq t+dt ) 1 m (B.6)
Note that the accumulated plastic strain increment ∆ε eq = ε eq t+dt -ε eq t is equal to Φdt. [START_REF] Lubarda | On the numerical algorithm for isotropickinematic hardening with the Armstrong-Frederick evolution of the back stress[END_REF]] or [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]], that

s t+dt -x t+dt σ eq t+dt = B pr |B pr | and ∆ε v = 3 2 ∆ε eq B pr |B pr | (B.9)
Based on relationship (6), the Huber von Mises stress can be related to the accumulated plastic strain increment ∆ε eq :

σ eq t+dt = (σ Y + k(ε eq t + ∆ε eq ) n ) 1 + ∆ε eq εo dt m . (B.10)
Therefore, combining Eqs (B.8) and (B.10), the accumulated plastic strain increment is shown to be solution of the non linear equation: F (∆ε eq ) = 0 where

F (∆ε eq ) = [(σ Y + k(ε eq t + ∆ε eq ) n ) 1 + ∆ε eq εo dt m + (3G + h)∆ε eq ] 2 -|B pr | 2 (B.11)
By a Newton-Raphson algorithm, the accumulated plastic strain increment ∆ε eq is evaluated. Consequently, ∆ε v is obtained from Eq. (B.9). Based on Eq. (B.4), s t+dt and x t+dt are updated.

It is also necessary to define the Jacobian matrix C J = ∂∆σ ∂∆ε for the implementation of the model in a UMAT subroutine. The main formula provided in [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]] are recalled. From ∆σ = C.∆ε -2G∆ε v , the Jacobian matrix is:

C J = C -2G ∂∆ε v ∂∆ε (B.12)
where C = (M e ) -1 is the fourth order tensor for linear elastic stiffness.

From Eq. (B.9), three derivatives need to be evaluated: ∂B pr ∂∆ε , ∂|B pr | ∂∆ε , ∂∆εeq ∂∆ε . With use of definition of B pr (B.5), it is easily shown that:

∂B pr ∂∆ε = 2GK ∂|B pr | ∂∆ε = 3G B pr |B pr | (B.13)
The term ∂∆ε eq ∂∆ε can be obtained without any difficulties from Eq. (B.11):

∂∆ε eq ∂∆ε = ∂(|B pr | 2 ) ∂∆ε H (B.14)
where H = ∂F ∂∆ε eq . One has to remember that the elastic predictor term |B pr | is related to the total strain increment, so that from Eq. (B.5):

∂(|B pr | 2 ) ∂∆ε = 12G 2 ∆e + 6G(s t -x t ) (B.15)
Combining Eqs. (B.14) and (B.15), ∂∆ε eq ∂∆ε is evaluated and the Jacobian matrix C J is obtained:

C J = ∂∆σ ∂∆ε = C - 6G 2 ∆ε eq |B pr | K + 9G 2 B pr |B pr | ⊗ B pr |B pr | ∆ε eq |B pr | -2 |B pr | H (B.16)
It is interesting to notice that from the work of [START_REF] Lubarda | On the numerical algorithm for isotropickinematic hardening with the Armstrong-Frederick evolution of the back stress[END_REF]] or [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]],

the update of the Cauchy stress tensor needs to solve only one non linear equation. Afterwards, the Jacobian Matrix can be evaluated in an analytical manner. The present Appendix is dedicated to the Perzyna type overstress model, but can be generalized to any other viscous flow law, see also [START_REF] Lubarda | On the numerical algorithm for isotropickinematic hardening with the Armstrong-Frederick evolution of the back stress[END_REF]] or [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]] for more information.

C Analytical solution for the Kroner MT model

Let us find the analytical solution for the Kroner MT model in the case of a two-phase composite subjected to proportional tensile loading. For such loading conditions, under assumed constitutive models, the stress in each phase σ, the overall stress Σ and strain rate per-phase ε are of the following form :

σ =               σ 11 0 0 0 σ 22 0 0 0 σ 22               Σ =               Σ 11 0 0 0 0 0 0 0 0               ε =               ε11 0 0 0 ε22 0 0 0 ε22               (C.1)
The overall strain rate Ė has the same structure as ε.

In the present paper, elasticity is linear isotropic. Linear kinematic hardening law is adopted with the same value of kinematic hardening modulus h for both phases. In addition, we restrict the development in this Appendix to the case without isotropic hardening. In that configuration, when plasticity develops, relation ( 6) for the inclusion and matrix phases takes the scalar form:

σ eq = 3 2 (s 11 -x 11 ) = σ Y 1 + εv 11 ε0 m , (C.2)
where s is the deviatoric part of the Cauchy stress tensor and, since the plasticity is volume preserving, the backstress x and the viscoplastic strain rate tensor εv are deviatoric by construction and take the form :

x =               x 11 0 0 0 -1 2 x 11 0 0 0 -1 2 x 11               , εv =               εv 11 0 0 0 -1 2 εv 11 0 0 0 -1 2 εv 11               (C.3)
The equation (C.2) can be rewritten as:

s 11 = x 11 + 2 3 σ Y 1 + εv 11 ε0 m (C.4)
Taking into account the linear evolution law (7) of the backstress, the rate of s 11 is specified as : . The elastic parameters of the phases are given in Table 1. The volume fraction of inclusions is c = 25%. As a reference, the slope for the purely elastic case is also shown.

ṡ11 = 2 3 h εv 11 + mσ Y

D Characterization of heterogeneity

This section depicts the way measures of heterogeneity are obtained at the inclusion level and at the finite element scale, see [START_REF] Czarnota | Modeling of the cyclic behavior of elastic-viscoplastic composites by the additive tangent Mori-Tanaka approach and validation by finite element calculations[END_REF]] for more information. For the sake of generality, relationships are provided in terms of a variable designated as κ and can be easily adapted to any field variables of interest (components of the strain, stress or backstress tensor). Note also that time dependence of the variable κ is intentionally omitted in this appendix.

The first scale of observation refers to κ I , corresponding to the mean value in each inclusion I of the finite element domain. Therefore, for a given time, κ I is expressed as:

κ I = N I p=1 κ p dv p N I p=1 dv p , (D.1)
where N I stands for the number of finite elements in one inclusion I, κ p is the corresponding local value inside the finite element p belonging to inclusion I and dv p the current volume of the element.

The second level of observation concerns the fluctuation of κ between all the finite elements belonging to the inclusion domain. By contrast with the first scale of observation, this second type of measure has no concern neither with the position of the finite element within the inclusion domain neither with the inclusion it belongs to. This measure relies on the identification, at each time t of the deformation process, of a probability density function Π κ (κ) defined on the interval [κ min , κ max ] by:

Π κ (κ) ≥ 0, ∀κ ∈ [κ min , κ max ] κmax κ min Π κ (κ) dκ = 1 , (D.2)
where κ min and κ max denote the minimum and maximum values of κ recorded in the FE calculation at time t. The probability of finding κ lying inside the interval [κ 1 , κ 2 ] is given by:

P κ (κ 1 ≤ κ ≤ κ 2 ) = κ 2 κ 1 Π κ (κ) dκ . (D.
3)

The two scales of heterogeneity are also characterized by the corresponding relative absolute deviations around the mean defined by:

r I κ = 1 N inc N inc I=1 |κ I -κ| |κ| , (

Fig. 1 .

 1 Fig. 1. Schematic representation of the composite. Spherical inclusions are dispersed in the matrix phase. The volume content of inclusion phase is c. The composite material is subjected to a macroscopic strain rate tensor Ė prescribed at the remote boundary. The inclusion problem for the Mori-Tanaka scheme is represented by an inclusion embedded in an infinite volume, made of matrix phase, subjected to the average matrix strain rate tensor εm .

Fig. 2 .

 2 Fig. 2. Presentation of the finite element model. A 3D representative volume element with 30 inclusions is adopted. The volume fraction of inclusions is c = 0.25. Periodic boundary conditions are adopted to prescribe the tension-compression cyclic loading.

Figure 3

 3 Figure3compares the predictions of the additive tangent MT model with finite element calculations based on a 30 inclusions RVE. The estimate provided by the MT scheme combined with the Kroner interaction law is also added. Note that by construction, only mean values can be evaluated with the proposed Mori-Tanaka scheme. From the overall response, it is clearly observed that the predictions of the additive tangent MT scheme are accurate, even from the first cycle. One can nevertheless notice little discrepancy during the first loading/unloading stage, see Fig.3 a). It must be also recognized that predictions of the Kroner model are correct too. At least for the overall behavior, the results obtained with the additive tangent interaction law is sligthly more precise, see Fig.3 a).

Fig. 3 .

 3 Fig. 3. Stress-strain response at a) the overall level, b) the inclusion level. The predictions of the present approach (Additive tangent MT) is compared to the one of the Kroner model (Kroner MT) and to the Finite Element calculations (FEM). The phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are displayed in Table1.

Fig. 5

 5 Fig. 5 displays time evolutions of the average stress, backstress and strain in each of the 30 inclusions.

  5 a) and b). In contrast, the heterogenity observed in the average longitudinal backstress (Fig.5 c-d), as well as the one related to the average strain (Fig.5 e-f) are well marked during the first loading/unloading stage.

Fig. 4 .

 4 Fig. 4. Time evolution of the a) stress and b) backstress for the matrix phase c) stress and d) backstress in the inclusion phase. The prediction of the proposed Mori-Tanaka scheme based on the additive tangent interaction law is compared to those of the Kroner model and of the FE results. The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table1.

Fig. 4 .

 4 Fig. 4.

Fig. 6 .

 6 Fig. 6. Heterogeneity between each individual inclusion and probability density functions at the finite element level of stress (a, b), backstress (c, d) and strain (e, f) at the end of the first tensile loading shown in Fig. 5 when t = 10s (or macroscopic strain component E 11 = 0.01). The volume fraction of inclusions is c = 0.25. The material parameters, provided in Table 1, are those of Fig. 4.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Values taken at each reversal time of relative absolute deviations around the mean, see Eqs (D.4-D.5) of Appendix D, and absolute difference (17) between the additive tangent MT and the FE mean values in the inclusion phase, related to a) stress b) backstress and c) strain. The volume fraction of inclusions is c = 0.25. The material parameters are those of Figure 6.

Fig. 9 .

 9 Fig. 9. Time evolution during cycles of the mean longitudinal stress (left axis) and backstress (right axis) in the a) matrix domain, b) inclusion phase. The two phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are those of Fig. 8.

Fig. 10 .

 10 Fig. 10. Stress-strain response a) of the RVE and b) of the inclusion phase. The prediction of the additive tangent Mori-Tanaka scheme is compared to the one obtained with a Kroner based scheme and to the FE results. The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1 except for the isotropic hardening parameter k which is set to k = 0GP a.

Fig. 11 .

 11 Fig. 11. Time evolution of the a) stress and b) backstress for the inclusion phase. The prediction of the additive tangent Mori-Tanaka scheme is compared to the one obtained with a Kroner based scheme and to the FE results.

Fig. 12 .

 12 Fig. 12. Stress-strain response for the overall, matrix and inclusion phases. The phases are elastic-viscoplastic. Isotropic and linear kinematic hardenings are disregarded. The volume fraction of inclusions is c = 0.25. The material parameters are displayed inTable 1 except for k and h which are set to zero.

Fig. 13 .

 13 Fig. 13. Macroscopic stress-strain response for three values of the kinematic hardening parameter h. The predictions of the proposed model are compared to the FE predictions. Only the first cycle is reported. The inclusion

FE

  model coincide. This provides an additional information concerning the plastic strain rate in the matrix. As already stated, because of linear kinematic hardening (7), a linear relation exists between the backstress rate and the plastic strain rate tensors. Therefore, from Fig.15 b), it is shown that the a)

  Fig. 15. a) Stress-strain response for the RVE and in the inclusion. b) Time evolution during cycles of the mean longitudinal stress σ m 11 (left axis) and backstress x m 11 (right axis) in the matrix. The matrix is elastic-viscoplastic with combined isotropic and linear kinematic hardening. Inclusions are elastic. The volume fraction of inclu-

Fig. 16 .

 16 Fig. 16. Stress-strain response at a) the macroscale, b) the inclusion level. The phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are displayed in Table2. The inclusion phase is soft.

Fig. 17 .Fig. 18 .Fig. 19 .

 171819 Fig. 17. Time evolution during cycles of the mean longitudinal stress (left axis) and backstress (right axis) in the a) matrix. b) inclusion. The phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are those of Fig. 16. The inclusion phase is soft.

  Fig.20again illustrates that the measures of the absolute deviation around the mean at the inclusion level and at the finite element level both follow the same trends for any of the observed field variables.

Fig. 20 .

 20 Fig. 20. Values taken at each point of reversal of relative absolute deviations around the mean defined in Eqs (D.4-D.5) and absolute difference (17) between the additive tangent MT and the FE mean values in the inclusion phase for the longitudinal a) stress b) backstress and c) strain. The volume fraction of inclusions is c = 0.25. The material parameters, provided in Table 2, are those of Fig. 16.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Stress-strain response at a) the overall level, b) the inclusion level. The phases are elastic-viscoplastic with linear kinematic hardening. No strain hardening is accounted for. The volume fraction of inclusions is c = 0.25. The material parameters are displayed in Table2except for k = 0GP a. The inclusion phase is soft.

  Fig. C.1. Evaluation of the normalized slope of the overall strain-stress curve in uniaxial tension after reaching the threshold σ Y in both phases as a function of kinematic hardening modulus h according to the Kroner MT model (relation (C.11)). The elastic parameters of the phases are given inTable 1. The volume fraction of

  D.4) for the inclusion level with κ I given by Eq. (D.1) and by: r κ = κmax κ min Π κ (κ)|κ -κ| dκ |κ| (D.5) for the finite element level. Note that κ stands for the average value of κ in the inclusion phase. So it corresponds respectively to the mean values defined in the present paper : σ i 11 , x i 11 and ε i 11 for the inclusion. N inc is the number of inclusions in the FE model (here N inc = 30). Time evolution jointly with close-ups of the heterogeneity of a-b) stress, c-d) backstress e-f) strain for all individual inclusions based on FE calculation. The average value in the inclusion phase obtained by FE and the prediction of the proposed Mori-Tanaka scheme based on the additive tangent interaction law are superimposed. The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1. 19 6 Heterogeneity between each individual inclusion and probability density functions at the finite element level of stress (a, b), backstress (c, d) and strain (e, f) at the end of the first tensile loading shown in Fig. 5 when t = 10s (or macroscopic strain component E 11 = 0.01). The volume fraction of inclusions is c = 0.25. The material parameters, provided in Table 1, are those of Fig. 4. 20 7 Values taken at each reversal time of relative absolute deviations around the mean, see Eqs (D.4-D.5) of Appendix D, and absolute difference (17) between the additive tangent MT and the FE mean values in the inclusion phase, related to a) stress b) backstress and c) strain. The volume fraction of inclusions is c = 0.25. The material parameters are those of Figure 6. 23 8 Stress-strain response at a) the overall level, b) the inclusion level. The phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are displayed in Table 1 except for the isotropic hardening coefficient involved in Eq. (5) whose value is k = 2.5GP a. 24 61 Time evolution during cycles of the mean longitudinal stress (left axis) and backstress (right axis) in the a) matrix domain, b) inclusion phase. The two phases are elastic-viscoplastic with combined isotropic and linear kinematic hardening. The volume fraction of inclusions is c = 0.25. The material parameters are those of Fig. 8. 25 10 Stress-strain response a) of the RVE and b) of the inclusion phase. The prediction of the additive tangent Mori-Tanaka scheme is compared to the one obtained with a Kroner based scheme and to the FE results. The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1 except for the isotropic hardening parameter k which is set to k = 0GP a. 27 11 Time evolution of the a) stress and b) backstress for the inclusion phase. The prediction of the additive tangent Mori-Tanaka scheme is compared to the one obtained with a Kroner based scheme and to the FE results. The volume fraction of inclusions is c = 0.25. The material parameters are provided in Table 1 except for the isotropic hardening parameter k which is set to k = 0GP a. 28 12 Stress-strain response for the overall, matrix and inclusion phases. The phases are elastic-viscoplastic. Isotropic and linear kinematic hardenings are disregarded. The volume fraction of inclusions is c = 0.25. The material parameters are displayed in

  Fig. 16. 37

  Table 1 except for k and h which are set to zero.

		1500
			MT h=10 GPa
			FEM h=10 GPa
	Σ 11 (MPa)	500 1000	MT h= 100 GPa FEM h= 100 GPa MT h= 1000 GPa FEM h= 1000 GPa
	Overall stress	-500 0
		-1000
				m=0.1, c=0.25
		-1500
		-0.01	-0.005	0	0.005	0.01
			Overall strain E 11

Table 2

 2

Table 1

 1 except for k and h which are set to zero. 29 13 Macroscopic stress-strain response for three values of the kinematic hardening parameter h. The predictions of the proposed model are compared to the FE predictions. Only the first cycle is reported. The inclusion is harder than the matrix, with material parameters given in Table1. The kinematic hardening modulus is

Acknowledgments

The research was partially supported by the project No. 2016/23/B/ST8/03418 of the National Science Centre, Poland, by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement Nº 777896 (Project QUANTIFY) and by the Agence Nationale de Recherche through the program Labcom LEMCI ANR-14-LAB7-0003-01.

A Time integration procedure for the Mori-Tanaka scheme under uniaxial tension/compression

The update of the mechanical state of the composite and the phases from time t to t + dt by the Mori-Tanaka scheme is carried out with a forward Euler scheme. More details can be found in [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and Mori-Tanaka schemes[END_REF]]. At time t, the Cauchy stress tensors in the matrix σ m (t) and in the inclusion σ i (t) are supposed to be known. The corresponding quantities for the backstress are also known quantities.

From the local behavior of the phases (see Eq. ( 8)), the viscoplastic strain rate tensors εv i (t) and εv m (t) are obtained. Thus, the fourth order tensor of viscoplastic compliance for both phases M v(tan) can be calculated, see Eq. ( 10). The corresponding elastic compliance tensors M e are given. It is interesting to notice that owing to the structure of the tangent additive interaction law, for the Mori-Tanaka scheme, only the viscoplastic M v(tan) m and elastic M e m compliance tensors are compulsory for the evaluation of the Eshelby tensor. This is an important advantage when compared to other theories of the literature where elastic-viscoplastic consistent tangent stiffness needs to be calculated.

From the interaction law (14) and from the incremental elastic law for all phases, the total strain rate in the inclusion phase is obtained at time t by Eq. ( 15). In the present work, the composite is under uniaxial tension/compression, at a prescribed longitudinal overall strain rate Ė11 . From the consistency equation, one gets a relationship for the strain rate component in the loading direction:

Owing to uniaxial loading tension/compression, the overall transversal stress components are Σ 22 = Σ 33 = 0 during all the loading, which is equivalent to Σ22 = Σ33 = 0. Since the macrostress is the volume average of the local stress in the inclusion and in the matrix, one gets from Eq. ( 16) :

Note that we have assumed in the model that all stress and strain rate tensors are diagonal ones, as a consequence of the isotropy of the phases.

From the Hooke's law (2), Eqs ( 1) and ( 15), the local strain rate tensors in the inclusion εi and in the matrix εm are found. From the knowledge of the local strain rate tensors, and the Hooke's law, the local stress tensors are updated. The backstress tensors are updated, via a forward Euler scheme too.

The process can be repeated for the next time increment.

B Numerical implementation

The numerical implementation for combined isotropic-kinematic hardening follows the classical radial return method. More details can be found for instance in [START_REF] Lubarda | On the numerical algorithm for isotropickinematic hardening with the Armstrong-Frederick evolution of the back stress[END_REF]] or [START_REF] Gomez | Damage mechanics constitutive model for Pb/Sn solder joints incorporating nonlinear kinematic hardening and rate dependent effects using a return mapping integration algorithm[END_REF]]. In the present case, linear isotropic elasticity is considered. We denotes by K the bulk modulus and by G the shear modulus.

Assume that the increment of the total strain tensor ∆ε is provided. One has to update the Cauchy stress tensor σ t+dt and the backstress x t+dt tensor, knowing their corresponding quantities at time t :

by assuming that εv 11 ≃ 0. Using relation (1) and the fact that εv is deviatoric, one obtains:

Note that (.) ′ denotes the deviatoric part of (.). From isotropic elasticity, ṡ11 = 2G εe ′ 11 , with G the elastic shear modulus of each phase. With Eq. (C.6), one gets :

Using the Kroner interaction law, for isotropic elastic behavior of the matrix phase, it is found that

After some standard manipulations of the Mori-Tanaka scheme, the link between the deviatoric parts of the overall strain rate tensor Ė and of the macrostress Σ is obtained:

where Ēh is the slope of the overall strain-stress curve when viscoplastic flow takes place. It is related to Ḡh and K in the same way as the elastic Young modulus E to the elastic bulk and shear moduli K and G. It has been verified that when h tends to infinity, the slope is approaching the one of the elastic regime. The material parameters are provided in Table 1. 18 60 C.1 Evaluation of the normalized slope of the overall strain-stress curve in uniaxial tension after reaching the threshold σ Y in both phases as a function of kinematic hardening modulus h according to the Kroner MT model (relation (C.11)). The elastic parameters of the phases are given in Table 1. The volume fraction of inclusions is c = 25%. As a reference, the slope for the purely elastic case is also shown. 52

List of Tables 1 Case of hard inclusions. Material parameters for the inclusion and matrix phases.

Viscoplasticity is described by the Perzyna-type law (6). Most of the parameters are taken from [START_REF] Pierard | Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions[END_REF]]. 14 2 Case of soft inclusions. Material parameters for the inclusion and matrix phases.

Viscoplasticity is described by the Perzyna-type law (6). 33