A double incremental variational procedure for elastoplastic composites with combined isotropic and linear kinematic hardening

We investigate the nonlinear behavior of elasto-plastic composites with isotropic and linear kinematic hardening. We first rely on the incremental variational principles introduced by Lahellec and Suquet [13]. We also take advantage of an alternative formulation, recently proposed by Agoras et al [1] for visco-plastic composites without hardening, which consists in a double application of the variational procedure of Ponte-Castañeda. We extend in this paper this approach to elastoplastic composites with combined linear kinematic and isotropic work-hardening. The first application of the variational procedure linearizes the local behavior, including hardening, and leads to a thermo-elastic Linear Comparison Composite (LCC) with a heterogeneous polarization field inside the phases. The second one deals with the heterogeneity of the polarization and results in a new thermo-elastic LCC with a per-phase homogeneous polarization field, which effective behavior can then be estimated by classical linear homogenization schemes. We develop and implement this new incremental variational procedure for composites comprised of linear elastic spherical particles isotropically distributed in an elastoplastic matrix. The predictions of the model are compared with results available in the literature for cyclic proportional and non-proportional loadings. New results for elasto-plastic composites with combined isotropic and kinematic hardening are also provided. They are in good agreement with the numerical computations we carried out, at both local and macroscopic scales. Preprint submitted to International Journal of Solids and Structures August 6, 2018


Introduction
The aim of this paper is to design a new incremental variational method able to simulate the behavior of composites made of elasto-plastic phases with both isotropic and linear kinematic hardening. Such composites exhibit reversible and irreversible effects when they are deformed. The presence of both these effects in the phases and their coupling generate some specific phenomena at the macroscopic level, such as the Bauschinger effect, which are still difficult to capture. Establishing the macroscopic behavior of such heterogeneous materials from properties and constitutive laws of the different phases is a classical but still open general problem in the mechanics of materials. Many methodologies and results have been proposed for linear-elastic composites, a revue of which can be found in Milton [19]. Many predictions for nonlinear composites have also been proposed and some renewal of these approches has been introduced in the late 80s with the introduction of variational formulations, e.g. Ponte Castañeda and Suquet [23]. Apart from the last ten years, most of the work published on this subject dealt with composites whose phases are governed by a single potential: the free-energy density in the case of nonlinear elastic materials or the dissipation potential in the case of purely viscoplastic materials. For such type of heterogeneous materials, the nonlinear homogenization theories rely on two steps: a linearization step followed by a homogenization one. The first step allows to approximate the nonlinear composite by a fictitious linear comparison composite LCC (Talbot and Willis [29], Willis [31], Ponte Castañeda [21,22], Suquet [27]) which is homogenized during the second step, by means of classical linear homogenization scheme. Earlier models made use of the first-order moments of the local fields only, i.e. their perphase averages, in the linearization schemes. Thanks to variational formulations and physical arguments, other authors highlighted the importance to incorporate also the fluctuations of the fields, e.g. by means of their second-order moments over the phases (Ponte Castañeda [21], Suquet [28], Ponte Castañeda and Suquet [23], Moulinec and Suquet [18]), or the discretization into subphases (Bilger et al. [2], Michel and Suquet [18]). These approaches significantly enhance the accuracy of the models. So far, the most sophisticated models consider both the first and second-order moments per phase to define the LCC (Ponte Castañeda [25]). The efficiency of such approaches has been supported by several theoretical and numerical studies (e.g. Idiart et al. [9], Idiart and Ponte Castañeda [10]). Some 10 years ago, a new direction of research has been opened by Lahellec and Suquet [13,14] who developed new incremental variational principles both at the local and global scales for composites with nonlinear hereditary behaviors. Similar variational principles had been previously developed by Mialon [17] for elastoplastic composites made of Generalized Standard Material (GSM) phases (Halphen and Nguyen [8]), and by Ortiz and Stainier [20] for elasto-viscoplastic materials in large deformations, but only at the local scale. The local variational principle of Lahellec and Suquet relies on the introduction of an unique potential, the condensed incremental potential, which accounts for both the conservative and dissipative effects of the composite behavior. This potential is built as the sum of the two potentials which define the local behavior in the framework of GSM, i.e. the free-energy density and the dissipation potential. From their local variational principle, Lahellec and Suquet derive an effective incremental variational principle for heterogeneous materials with hereditary behavior thanks to which it becomes possible to extend to composites with hereditary behavior most of the ideas developed in the variational framework of nonlinear homogenization for composites governed by a single potential. Following this way, Lahellec and Suquet [13] applied the key idea of the variational procedure introduced by Ponte Castañeda [21] to their new incremental variational principles to estimate the behavior of nonlinear composites made of elasto-viscoplastic phases. For that, they made use of two approximations. The phase behavior is linearized by introducing a uniform secant viscosity per phase through the variational procedure. Furthermore, the heterogeneous plastic strains at time t n are approximated by an effective homogeneous internal variable in each individual phase. These two approximations lead to the definition of a thermo-elastic LCC (Ponte Castañeda [21]) whose effective behavior can then be derived by any appropriate linear homogenization scheme appropriate to the considered microstructure. This procedure, called the Effective Internal Variables (EIV) approach, was first applied to nonlinear elastoviscoplastic composites with neither threshold nor hardening and leads to accurate estimates of the macroscopic behavior. In 2013, Lahellec and Suquet [16] introduced modified incremental variational principles in a rate form in order to get estimates of the local and global behavior for composites made of elasto-viscoplastic phases with both isotropic and linear kinematic hardening. Their new approach, called the Rate Variational Procedure (RVP), still relies on two approximations. The first one consists in using the variational procedure to linearize the constitutive laws and provides again a secant approximation of the behavior. For the second one, at the difference with Lahellec and Suquet [13], it is now the heterogeneous stress field -and not the heterogeneous plastic strain field -at time t n which is approximated by effective homogeneous internal variables in each individual phase. For that, they make use of a similar method than the one proposed by Lahellec et al. [15], which is also based on the variational procedure but relies on a new definition of a LCC with uniform coefficients per phase which depend on the first and second-order moments of the stress field. In order to assess their new model, Lahellec and Suquet [16] carry out Fast Fourier Transform (FFT) simulations on a heterogeneous material made of an elasto-(visco)plastic matrix reinforced by 50 elastic isotropic spherical particles randomly distributed. A good agreement is observed between the model predictions and the FFT simulations for elastic ideally-plastic matrix or elasto-plastic matrix with isotropic or linear kinematic hardening under radial and non-radial loadings. Similarly, Boudet et al. [3], relying on the initial variational principles developed by Lahellec and Suquet [13], extended their EIV formulation to elasto-(visco)plastic composites with both threshold and local isotropic and linear kinematic hardening. The results given by their approach are very close to the numerical simulations of Lahellec and Suquet [16]. It should be noted that Brassart et al. [5,6] also developed an incremental variational procedure relying on the variational principle introduced by Ortiz and Stainier [20] and obtained accurate predictions of the effective behavior of elasto-(visco)plastic composites with local isotropic hardening. More recently, Wu et al. [32,33] developed a new incremental homogenization scheme for elasto-(visco)plastic composites which, instead of relying on variational principles, is based on incremental secant tensors, homogeneous by phase, and evaluated by taking into account the residual stress and strain fields. In their original version [32], the linearization is performed by introducing only first-order moments of the local fields while in their second version use is made of second statistical moments. Lastly, Agoras et al. [1] proposed an alternative formulation to the RVP model of Lahellec and Suquet [16] for elasto-viscoplastic composites without hardening. Their approach consists in a double application of the variational procedure to the incremental condensed stress potential. The first application linearizes the local behavior and leads to a LCC with non-uniform phase properties, namely a heterogeneous eigenstrain within the phases. The second application, which deals with this heterogeneity of phase properties, makes use of the method developed by Lahellec et al. [15] to approximate a LCC with non-uniform phase properties by a LCC with homogeneous properties by phase. Note that, roughly speaking, the main difference between the new approach proposed by Agoras et al. [1] and the Lahellec and Suquet's ones [13,16] lies in the fact that the two main difficul-ties of the homogenization problem to be solved, that is the linearization of the local behavior and the accounting for the heterogeneity of the LCC introduced by the local fields at time t n (the plastic strain field [13] or the stress field [16]), are dealt with sequentially by Agoras et al. [1] instead of being addressed simultaneously as [13,16]. Such a treatment allows to define more clearly the main two problems to be solved and to handle them in a more efficient and systematic way. Agoras et al. [1] apply their approach to particulate composites made of elastic ideally-plastic polymer or metal matrix reinforced (or softened) by spheroidal elastic particles. For the ideally-plastic metal matrix reinforced by spherical elastic particles, Agoras et al.'s model and the Lahellec and Suquet's FFT simulations are in good agreement even if the stress averages over the inclusions, as well as the stress-fluctuations within the matrix, are significantly overestimated, as it is also the case for both the RVP formulation of Lahellec and Suquet and the model of Boudet et al. [3].
The approach presented in this paper is based on the key idea presented by Agoras et al. [1] to handle sequentially the linearization of the local behavior and the accounting for the heterogeneity of the LCC. We propose to extend this idea, initially applied to elasto-viscoplastic composites without hardening, to elastoplastic composites with isotropic and linear kinematic hardening. As Agoras et al. [1], we still make use of a double application of the variational procedure to the condensed local incremental potential. However, instead of using a stressbased energy as in [1], our derivation is formulated in the strain-based free-energy framework. The first application of the variational procedure aims at linearizing the local behavior, including hardening, and leads to a thermo-elastic LCC with heterogeneous polarization field inside the phases. The second one deals with the heterogeneity of the polarization by making use of the method developed by Lahellec et al. [15]. It results in a thermo-elastic LCC with homogeneous polarizations inside the phases, which is estimated by classical homogenization schemes. The resulting model which, for the sake of simplicity, will be referred to in the sequel as the Double Incremental Variational (DIV) procedure, is then compared with the RVP approach and FFT simulations for several composites under radial and non-radial loadings [16]. For the purpose of such comparisons, the composites under consideration are made of either an elastic ideally-plastic matrix or an elasto-plastic matrix with isotropic hardening or linear kinematic hardening reinforced with spherical elastic particles. In this context, the DIV formulation is also compared to other results of the literature such as Brassart et al. [6] derived for elasto-plastic composites with isotropic hardening under axial tensile loading. Lastly, new results for elasto-plastic composites with combined isotropic and lin-ear kinematic hardening are also provided and show that the new formulation is in good agreement with Finite Element (FE) simulations.
The structure of the paper is as follows. Inspired by the work of Lahellec and Suquet [13], Section 2 briefly recalls the local and global variational principles. It also defines the elasto-plastic behavior of the phases which is described by the J 2 theory of plasticity combined with both linear kinematic harding and nonlinear isotropic hardening. Section 3 presents the development of the DIV procedure based on the work of Agoras et al. [1] and is divided into four subsections. The first one explains how to linearize the behavior by means of the variational procedure and to derive a LCC with a heterogeneous polarization within phases while the second one presents how to approximate this heterogeneous LCC by a "classical" one with a homogeneous polarization per phase. In the third subsection, the DIV formulation is then particularized to the case of isotropic elasto-plastic composites reinforced by spherical linear elastic particles randomly distributed within the matrix. The algorithm implementation is described in the last subsection as well as the difficulties met during the control of the convergence (numerical accuracy and non-uniqueness of the solution). Finally, Section 4 deals with the comparison between the DIV approach and the data found in the literature as well as new FE simulations.
The tensor notation used herein is a fairly standard one. The order of the tensors are clear when taken in the context. Products containing dots denote summation over repeated indices. For example, L : ε = L ijkl ε kl e i ⊗ e j and E :: F = E ijkl F klij where e i (i = 1, 2, 3) is a time-independent orthogonal Cartesian basis and the operation ⊗ denotes the classical tensorial dyadic product.

Local behavior and incremental variational principles
This section is inspired by the work made by Lahellec and Suquet [13]. In what follows, we summarize the main steps and equations describing the local behavior and the local and global variational principles.

Local behavior
We consider a Representative Volume Element (RVE) Ω of a composite made of N phases occupying domains Ω (r) (r = 1, . . ., N ) such that Ω = ∪ N r=1 Ω (r) . The volume fraction of phase r is denoted c (r) with c (r) = |Ω (r) |/|Ω|. The phase 6 distribution is characterized by the characteristic functions χ (r) such as The phases are generalized standard materials (GSM) having an elasto-plastic local behavior described by the conventional J 2 theory of plasticity combined with both linear kinematic hardening and nonlinear isotropic hardening. This corresponds to a material with internal variables α = (ε p , p) describing the irreversible phenomena, where ε p is the plastic strain and p the cumulated plastic strain, and with two convex potentials. The first one is the free-energy density w r (ε, ε p , p) where ε denotes the local strain with L (r) the elasticity tensor, H (r) the kinematic hardening fourth-order tensor andŵ r (p) a scalar function which depends on the cumulated plastic strain p according toŵ where R(p) is a scalar function characterizing the isotropic hardening. The stress σ and the thermodynamic forces A α = (A ε p , A p ) associated with the internal variables α = (ε p , p) are given by the state laws: with X the so-called back-stress.
The second potential is the dissipation potential ϕ (r) (α) given by (e.g. [11], [7]): with and which is the indicator function of the convex set C defined as The evolution equations are then given by Remark 1. The Legendre-Fenchel transform gives access to the dual potential ϕ * (r) (A ε p , A p ) of the dissipation potential ϕ (r) which is readily shown to read such that ϕ * (r) is the indicator function Φ C * of the convex set if one chooses to define R(p) such that R(0) = 0.,

Local and global variational principles
As Lahellec and Suquet [13] we work with a "total" formulation, i.e. depending on ε, ε p and p, instead of a "rate" one, that would depend onε,ε p andṗ as in Lahellec and Suquet [16].

Local variational principle
From (4) and (9) we classically obtain the constitutive equations for GSM: The considered time interval [0, T ] is discretized into N time steps, not necessarily of identical duration and characterized by the set {t 0 = 0, t 1 , ..., t n , t n+1 , ..., t N = T }. For convenience, the value f (t n ) of a function f at time t n is denoted f n . At time t = t n we define the time increment ∆t = t n+1 − t n . As explained in Lahellec and Suquet [13] the time derivative in (12) can be approximated by the following equations using an implicit Euler-scheme In the sequel, to simplify the notations, we will omit the index n + 1 for the variables computed at time t n+1 (i.e. ε = ε n+1 ).
As [13] we define the local condensed incremental potential w ∆ and the incremental potential J as Let us emphasize on the fact that the dependance on x of J (r) results from ε p n and p n which are respectively a second-order tensor and a scalar fields which result from the previous time step. The second equation in (13) is the Euler-Lagrange equation of the infimum condition in (14). In addition, as consequence of this condition, the stress at any points x ∈ Ω is given by The local problem which determines the local stress and strain fields in the RVE as a function of the prescribed macroscopic loading conditions reads then where the fields σ, ε and α depend on x and t. The symbol . (resp. . (r) ) defines the spatial average over Ω (resp. over Ω (r) ), E the macroscopic strain and ∂Ω denotes the boundary of Ω. Assuming macrohomogeneous boundary conditions, the details of the later do not need to be specified.

Macroscopic variational principle
A variational representation of the local problem can be obtained by noticing that ε is solution of the following minimum potential energy principle Defining the effective condensed incremental potential bỹ Lahellec and Suquet [13] show, thanks to the relation σ = ∂w ∆ /∂ε and Hill's lemma, that the macroscopic stress Σ = σ is given by 3. Double application of the variational procedure to elastoplastic composites with hardening The incremental variational procedure, proposed by Lahellec and Suquet [13], leads to the definition of a linear comparison composite (LCC) with homogeneous polarization per phase, characterized by a free energy w 0 . The variational procedure is applied once to deal with both the nonlinearity of the phases and the heterogeneity of the local incremental potential J within the phases -this heterogeneity being due to the dependance of J on α n (x) = (ε p n (x), p n (x)). In this method, when defining the linearized incremental potential J 0 , it is not easy to exhibit the adequate closed-form expression of J 0 , since two stages, namely the linearization of the behavior and the handling of the heterogeneity of J, are melted. Agoras et al. [1] introduced a more systematic method in which the linearization of the behavior and the heterogeneity of the incremental potential within the phases are addressed sequentially in two separate steps. The first step makes use of the variational procedure of Ponte Castañeda [21] to obtain a LCC with heterogeneous eigenstrains within the phases. The second one which also relies on the variational procedure makes use of the method proposed by Lahellec et al. [15] to reduce the resulting problem to a different LCC with now homogeneous properties.
In the present study, we reformulate the key idea proposed by Agoras et al. [1], which consists in addressing sequentially, first the linearization of the local behavior and then the heterogeneity of the resulting LCC, and which was initially applied to ideally-plastic phases, to the context of elastoplastic composites with hardening phases. Both isotropic and linear kinematic hardening are considered. To clearly distinguish our approach from previous works of the literature adressing the same issue, we first emphasize that our formulation relies on the initial variational incremental principle introduced by Lahellec and Suquet [13] in a total form and not on the modified version proposed by Agoras et al. [1] in a rate form, i.e. based on the strain rateε instead of the total strain ε 1 . We also stress out that our approach is established in a primal form, i.e. based on w(ε) and ϕ(α), while a dual formulation was considered by Agoras et al. [1]. Lastly, the way to deal with the isotropic and linear kinematic hardening is different from the one considered in Lahellec and Suquet [16] as it will be shown below.

Linearization of local behavior
Owing to the non quadratic character of the dissipation potential ϕ (r) the incremental potential J defined by (14) is difficult to homogenize. To bypass this difficulty a linearized incremental potential J (r) L is introduced in order to approach the dissipation potential ϕ (r) by a quadratic function ofε p which depends on a viscosity η (r) ε p uniform in phase r: A key idea of the variational procedure is to add and subtract to the potential J the potential J L , i.e. J = J L + J − J L , such that the first term J L can be homogenized by classical methods for linear media thanks to the quadratic part, while the difference J − J L can still be evaluated separately.
3.1.1. General form of the proposed estimate of the effective incremental condensed potential To begin with, let us introduce which can be rewritten, thanks to (5), as By injecting (21) and (22) into the variational formulation (18), using the definition (14) and noticing that the minimization over ε p and p can be performed separately so that the infimum can be taken out from the average, one gets: In (22), Φ C (ε p ,ṗ) = +∞ when g(ε p ,ṗ) > 0. As a consequence, the infimum in (23) is obtained under the condition g ε p −ε p n ∆t , p−pn ∆t ≤ 0. It is convenient to introduce ∆J bis as: such that (25) As in Lahellec and Suquet [13], rigorous upper bound ofw ∆ can be obtained by taking the supremum of ∆J bis over (ε p , p) under the constraint g ≤ 0: The constraint g(ε p ,ṗ) ≤ 0 in (25) with the definition (8) can be readily replaced by Moreover, from (24) ∆J bis can also be seen as a function of (ε p ,ṗ). It follows that Eq. (26) takes the form ∆J bis (.,ε p ,ṗ, {η ε p }) .
(28) However, as pointed out by Lahellec and Suquet [12,13], this upper bound can be too stiff in certain cases. Ponte Castañeda and Willis [24] and Ponte Castañeda [25] have observed for similar situations that this upper bound can be replaced by a sharper estimate through relaxing the supremum condition by a stationarity one. Doing so, one gets the following estimate forw ∆ : 13 Then, a last optimization with respect to (w.r.t.) the set {η ε p } leads to the final DIV estimate forw ∆ defined as

Stationary conditions
In this section, we develop all the stationary conditions which appear in expressions (29) and (30). First, we compute the stationarity of ∆J bis w.r.t. (ε p ,ṗ), thus providing the definition of the viscous modulus η (r) ε p (see Eq. (33)). Then, we make use of the stationarity ofw ∆ w.r.t. η , the former denoting the value ofε p at the stationarity of ∆J bis and the later the second-order moment ofε p in phase r. Finally, stationarity of J L w.r.t. (ε p , p) yields a closed-form expression of ε p as function of ε (see Eq. (42)) and provides an interpretation of the viscous modulus η Stationarity of ∆J bis .
The stationarity of ∆J bis at point x w.r.t.ε p andṗ is determined by applying the Karush-Kuhn-Tucker optimality conditions (KKT) where ∂ denotes the variation inε p andṗ. Making use of definitions (8) and (24) , for each x in phase r one gets In this system λ,ε p eq andṗ are unknown fields which may depend on the position x in phase r. However, φ (r) and R (r) take the same form in each point x in phase r; As a consequence of our modeling options, η (r) ε p is also uniform in phase r. If we assume in addition that p n is uniform in phase r, then the solution of system (32) is uniform per phase. Finally, since p 0 at the begining of our loading history is uniformly zero, it is shown by recursion, using the relation p = p n + ∆tṗ, that all these fields are uniform in phase r. In particular,ε p eq will be denoted ε p (r) ; likewise p (r) ≡ p, p (r) n ≡ p n and λ (r) ≡ λ. We now consider the two following cases for the lagrange multiplier: λ (r) > 0 and λ (r) = 0.
Combining with Eq. (32) 1 and Eq. (32) 2 leads to the following secant-like formula Case 2 : a second possibility is λ (r) = 0. According to Eq. (32) 2 , this would lead to R (r) (p) = 0 which appears to be non physical when isotropic hardening occurs. Accordingly, this case will not be considered.
Owing to Eq. (33), we are now able to calculate ∆J bis under the stationarity conditions. Indeed, noticing by means of Eq. (33) that η is an inversible function of ε p (r) since R (r) (p) is an increasing function, and denoting ε p (r) η This expression will be helpful in section 3.1.3.

Stationarity ofw
Making use of Eqs. (20), (24) and (29), the stationarity condition ofw V AR (35) Due to Eq. (33) the last term of (35) vanishes. It follows where a (r) = ι a d : a d (r) denotes the second-order moment of a deviatoric part second-order tensor field a with ι = 2 3 for a strain field and ι = 3 2 for a stress field. Accordingly, the optimal viscosity η (r) and p (r) = ∆tε p (r) (38) It is a sum of two viscosities respectively associated with perfectly plastic behavior for η (r) φ,sct and isotropic hardening for η The constraint h(ε p , p) ≤ 0 for h defined by Eq. (27) has to be satisfied for all x ∈ Ω (r) since ε p n and p n are heterogeneous fields. To circumvent this difficulty, this constraint is relaxed into N per-phase average conditions only, i.e. h (r) ≤ 0 for r ∈ {1, ..., N }. Stationarity of J L w.r.t. (ε p , p) under these new constraints is then looked for by expliciting again the KKT optimality conditions (see AppendixA) where λ (r) are N scalars and ∂ denotes the variation w.r.t. ε p and p. Making use of (20) and (27) together with (39), one gets where K = I − J denotes the deviatoric fourth-order isotropic projector, I the identity fourth-order symmetric tensor and J the spherical fourth-order isotropic projector whose components are respectively I ijkl = 1 2 (δ ik δ jl + δ il δ jk ) and J ijkl = To close the set of equations of the model, we now have to find a relation between ε p and ε. To this end, two cases are considered for the Lagrange multiplier: λ (r) > 0 and λ (r) = 0.
= 0 by means of Eq. (40) 3 . Therefore ε p = ε p n for each x. Such a situation corresponds to the elastic regime. It will be fully handled later in Section 3.3.2.
Case 2: λ (r) = 0. This situation corresponds to the plastic regime since (40) 2 does not constraint p − p n to be zero, as in the later case. Eq. (40) 1 with λ (r) = 0 reads from which one gets the expected expression for the plastic strain field, namely Others information can be obtained from Eq. (41). Indeed, owing to Eq. (4) 2 , Eq.
(41) reads also where σ d denotes the deviatoric part of the second-order tensor σ. As shown by Eqs. (37), (38) and (43), η (r) ε p can be interpreted as the modified secant viscous modulus associated with the elastoplastic behavior with isotropic hardening and computed at the second-order moment of the plastic strain rateε p (r) . From Eqs. (43) and (37), it is readily seen that Thus, for the DIV approach the plastic yield criterion (11) associated with the local behavior is satisfied when applied to the second-order moment of the local fields. Indeed, the local stress (σ − X) eq is approximated by σ − X (r) while the hardening stress R (r) (p) is evaluated at p (r) whose evolution is provided by the second-order moment of the plastic strain rateε p (r) (see Eq. (38)).

Determination of the LCC with heterogeneous intraphase polarization
To proceed further, ε p is replaced by its expression (42) in the definition (20) of J (r) L , thus allowing to define a thermoelastic LCC of per-phase free-energy w (r) L given by with (46) Such an energy corresponds to a LCC with heterogeneous intraphase fields within phase r, the polarization τ  L (x), depending on x through ε p n (x). Accordingly, the effective energy of this LCC is defined bỹ From Eqs. (29), (34) and (36), it follows 2 with f (r) defined by Eq. (34) 2 . Eventually, taking into account the stationarity conditions on η (r) wherew L (E) denotes the value ofw L (E, {η ε p }) calculated for the optimal set {η ε p }.

3.2.
Homogenization of the LCC with intraphase heterogeneous polarization 3.2.1. General procedure As seen above, the LCC of energy w L corresponds to a thermoelastic composite with heterogeneous intraphase polarizations. Although its effective be-haviour could be derived by costly full-field computations, it can not be obtained by classical linear mean-field homogenization schemes since they only apply to thermoelastic composites with homogeneous intraphase polarization. One way to homogenize the thermoelastic LCC of energy w L is to resort to the TFA method of Dvorak [4] which, as a first approximation, replaces the heterogenous polarization τ (r) L (x) with a per-phase constant field. However, it is well known that the TFA method provides much too stiff estimates of the effective behavior, e.g. [18]. In order to deal with the heterogeneity of the polarizations, Lahellec et al. [15] proposed a new and efficient method based on the variational procedure of Ponte Castañeda [22] which approximates the energy w L by an energy w 0 associated with a classical LCC with homogeneous intraphase polarizations τ L such that it becomes possible to derive the effective behavior as well as the statistics of the local fields by means of classical linear homogenization schemes 3 . This idea was first implemented by Agoras et al [1] for viscoplastic composites without hardening. In what follows, we implement it for elasto-plastic composites with isotropic and linear kinematic hardening.
Applying the key idea of the variational procedure of Ponte Castañeda [22], the energy w (r) L is rewritten as where w 0 corresponds to the energy of a classical thermoelastic LCC with homogeneous intraphase polarizations defined by (51) Following the procedure of Lahellec et al. [15] which is briefly recalled in AppendixB, on can obtain an estimate of the effective energyw L (E) given by the following 3 In the initial work of Lahellec et al. [15], the problem to be solved was more general since not only τ .
(52) withw and for which the notation ∆a (r) (.) = a (r) A straightforward calculation shows that the stationarity condition w.r.t. f where ε 0 corresponds to the strain field within the LCC defined by w Remark 2. As shown by Eq. (46) the tensor τ (r) L is purely deviatoric. Moreover, the nonlinear behavior being only sensitive to the deviatoric part, we assume a purely deviatoric polarization τ (r) 0 . Accordingly, ∆τ (r) is also deviatoric as well as ∆L (r) : ε 0 .
Eq. (55) 1 amounts to 21 where a (r) =< a > (r) denotes the average of any field a(x) over the phase r.
Note that τ L (r) can be easily obtained from ε p n (r) since, as shown by Eq. (46), the polarization τ (r) L linearly depends on ε p n . Making use of (55) 1 , Eq. (55) 2 can be rewritten as where C (r) (s) is the covariance tensor of the intraphase fluctuations of a secondorder tensor field s within phase r, defined by So Eq. (57) becomes where the second equality in Eq. (59) is obtained by using the property C (r) τ (r) resulting from the uniformity of τ The effective behavior of the nonlinear composite can now be computed from the estimate (52) as In Eq. (60), the approximation comes from Eq. (49), the third equality from the stationarity ofw L w.r.t. L Lastly, the first and second-order moments of the local fields of the nonlinear composite are respectively approximated by the ones of the LCC with heterogenous intraphase polarization defined by the energy w L . Moreover, Lahellec et al. [15] shown that the first and second-order moments of the strain field within the LCC w (r) L can be estimated by those of the LCC w (r) 0 . So, one gets where ε, ε L respectively denote the strain fields of the nonlinear composite and of the LCC of energy w L . The same approximations are used for the first and second-order moments of the stress field.

Elasto-plastic composites with isotropic phases
In this study, the behavior of the phases described in Section 2.1 is also assumed isotropic. Accordingly, the elastic tensors of moduli and the linear kinematic hardening tensor H (r) read Due to the isotropic behavior of the phase and for simplicity, the elastic tensors L As shown in Remark 2, ∆L (r) : ε 0 is deviatoric thus implying k (r) Under these conditions, (46) rewrites as Similarly, Eq. (42) becomes where ε d denotes the deviatoric part of ε. To fully determine L should be determined (see (37)). Recalling thatε p = ε p −ε p n ∆t , it readily follows from (66) thaṫ (67) Equation (67) requires to compute ε d : ε p n (r) . To this end, we introduce (62) 2 into (43) such that : from which the second-order moment of σ d can be derived. Introducingε p = ε p −ε p n ∆t together with Eq. (66) into Eq. (68) allows to isolate ε d : ε p n (r) . One gets .
(69) As seen above, the second-order moments ε (r) , σ (r) are approximated by ε 0 (r) , σ 0 (r) , the latter being derived respectively by (D.3) 2 and the following expression which is deduced from 3.3. Application to an elasto-plastic matrix reinforced by linear elastic particles We now specialize the DIV approach to the case of two-phase reinforced elasto-plastic composites, made of spherical linear elastic particles randomly and isotropically distributed inside an elasto-plastic matrix. Further, the effective properties of the LCC as well as the field statistics are evaluated by the Hashin and Shtrikman lower bound (see AppendixD), which is known to provide accurate results for such type of microstructure in the linear case when the volume fraction is moderate. In what follows, subscripts (1) and (2) denote the inclusion and matrix phase, respectively.
The first-order moments of the strain ε 0 over the phase read as where the localization tensors A : where P 0 corresponds to the classical Hill's tensor associated with the matrix phase of the homogeneous LCC for an isotropic distribution of spherical inclusions. To determine the first-order moments ε 0 (r) , the tensors L have to be computed first. These latter depend both on the first and second-order moments of the plastic strain field ε p n at time t n (see Eqs. (59) and (65)), which are already known from the previous step, and also on the second-order moments of the plastic strain rateε p (r) which remain to be determined. Once the ε 0 (r) are known, the effective stress Σ defined in (60) is given by When using the Hashin Shtrikman lower bound to evaluate the strain and stress fields of a thermoelastic linear composite, it is well known that the latter are uniform inside the inclusions such that C (1) (ε 0 ) = 0. Therefore, Eq. (59) for r = 1 amounts to Nevertheless, nothing is known about L L . This tensor is therefore indeterminate. To raise this indeterminacy we adopt the assumption used in the works of Lahellec et al. [15] and Agoras et al. [1] L (1) so that Eq. (56) rewrites as 3.3.1. Summary of the nonlinear system In this section, we summarize the equations to be satisfied by the DIV approach. Let us first recall that, as shown by Eq (64), the spherical part of the local behavior in the LCCs defined by w L and w 0 is indentical to that of the nonlinear composite.
In the inclusion phase, sinceε p In the matrix phase, the expressions of µ 0 , µ L and τ L which characterize the LCCs of free energy w 0 and w L are obtained as follows. The secant modulus µ .
Likewise, applying Eq. (59) to isotropic LCC leads to the following expression for µ : Note that once η ε p is known -or equivalentlyε p , see Eqs. (37) -the numerator in the square root of Eq. (80) is also known since τ L depends linearly on ε p n whose first and second-order moments are provided from the previous step. On the other hand, the first-order moment ε d 0 (2) can be derived by Eqs. (72) and (73) Of course, for isotropic composites both these moments depend on µ 0 and η (2) ε p . Note that the dependence of the problem w.r.t. the last variable η (2) ε p appears in expressions (65) of µ (2) L and τ (2) L . As seen above, the determination of η (2) ε p by Eq. (37) amounts to computeε p (2) . For that, use is made of Eq. (67) which itself requires to compute the second-order moment ε d : ε p n (2) . The latter is obtained from Eqs. (69)   L shows that τ (2) 0 also depends only on µ (2) 0 and η (2) ε p . Therefore, the nonlinear problem to be solved and defined by Eqs. (37), (79) and (80) now reduces to a system of two nonlinear scalar equations with unknown µ (2) 0 and η (2) ε p : After resolution of this system, the remaining parameters µ L , τ  79) and (67), respectively.

Case of the elastic regime
In the elastic regime, ε p = ε p n andε p = 0. The second equation of system (82) implies that η (2) ε p tends to +∞, so that this system reduces to a single nonlinear scalar equation F 1 µ (2) 0 , +∞ = 0. To find the required mechanical quantities during such an elastic loading and unloading, their limits when η (2) ε p tends to +∞ have to be computed. Doing so, we obtain in particular µ (2) L = µ (2) and τ (2) L = −2µ (2) ε p n . Furthermore, as we saw in Section 3.1.2, the plastic threshold is reached when Eq. (44) is satisfied in the matrix phase. Accordingly, to check whether we are in the elastic regime, the second-order moment σ − X (2) should satisfy the following inequality which is adopted as the yield surface for the matrix in our approach. For that, we need to compute the second-order moment σ − X (2) . As seen in Section 3.2.1, the first and second-order moments of the local fields in the nonlinear composite and in the LCC of energy w L are evaluated by their counterparts in the homogeneous LCC, e.g. see Eq. (61), so that where σ (2) is given by (70) and X = a (2) ε p n (2) because of Eq. (62) 2 and ε p = ε p n . Moreover, since σ d = 2µ (2) where the second equality comes from the expression (69) of ε d : ε p n (2) when η (2) ε p tends to +∞. Finally, σ − X (2) reads as

Numerical implementation of the model
We now consider the numerical implementation of the DIV model. First, we describe the structure of the proposed algorithm and then discuss its convergence to the appropriate solution.

Presentation of the algorithm
The aim of the algorithm is to determine, for a given loading history defined in terms of prescribed macroscopic strain E n at all times t n (n = {0, ..., N }), the effective stress Σ n , the first and second-order moments of the local fields σ n , ε n and ε p n in the phases as well as the value of p n . This problem will be solved iteratively by determining the solution at t n+1 from the known solution at time t n . As already noted, the index n + 1 will be obmited on the quantities defined at time t n+1 . This algorithm can be applied to an elastoplastic matrix which may be elastic ideally-plastic, with linear kinematic hardening, isotropic hardening or both isotropic and linear kinematic hardenings. For all these cases, the algorithm is the same and its structure goes as follows: 1. Elastic prediction: The matrix is assumed to follow an elastic regime such that η (2) ε p tends to +∞ according to system (82) which reduces to the nonlinear scalar equation F 1 µ  0 of the LCC with homogeneous intraphase polarization. Lastly, the effective free-energy as well as the statistics of the local fields -that is the first and second-order moments of ε and σ -are computed (see AppendixD) 4 . The effective stress is finally obtained by Eq. (60).
2. Check of the yield criterion: Now that all the coefficients of the LCC of energy w 0 as well as the statistics of the local fields have been determined, it becomes possible to evaluate σ − X (2) by means of Eq. (86) and to check if condition (83) is satisfied with a strict inequality. If it is the case, the matrix is well in the elastic regime and we can go to the next time step. If not, the matrix follows a plastic regime, discussed now.
3. Plastic regime : Unlike the elastic regime, η (2) ε p takes now finite value. The nonlinear problem to be solved is described by system (82) whose unknowns are µ (2) 0 and η (2) ε p . The procedure fsolve from Matlab is used again to find the two roots of this system. Once the solution of the nonlinear system (82) is obtained, the heterogeneous and homogeneous LCCs, as well as the effective free energyw 0 are determined as explained in Section 3.3.1. It is then possible to calculate the first and second-order moments of ε, σ and ε p , the cumulative plastic strain p (2) as well as the effective stress Σ.

Numerical accuracy
Let us now consider the influence of the time step discretization on the model predictions. To this end, simulations have been carried out on the same composite material as studied by Lahellec and Suquet [16] for different values of the time step ∆t. This composite is made of an ideally-plastic matrix reinforced with spherical elastic particles isotropically distributed. The considered material parameters are given in Section 4.1.2 (see (93)). The composite is submitted to a proportional isochoric macroscopic strain E(t) = E 33 (t) − 1 2 (e 1 ⊗ e 1 + e 2 ⊗ e 2 ) + e 3 ⊗ e 3 . Figure 1a reports the evolution of the axial macroscopic stress as a function of the axial macroscopic strain at four different time steps. Figures 1b and Figure 1c provide zooms on the transition zone between elastic and plastic regimes and on the fully plastic regime. As seen on these Figures, the four curves corresponding to the different time steps are almost identical, except for the highest time step increment for which a slight difference is observed in the transition zone.
For completeness, the evolutions of the stress fluctuations within the matrix, classically quantified by C (2) (σ) :: K -the stress covariance tensor C (2) (σ) being defined by Eq. (58) -are also depicted on Figure 1d for the different time steps. Again, the corresponding four curves are almost identical. Accordingly, the time step for all the simulations is fixed at ∆t = 0.1s.

Computational issues
As seen in Section 3.4.1, the determination of the effective behavior associated with the DIV approach as well as the statistics of the local fields requires to solve the nonlinear system (82) with two unknowns. Due to the form of Eq. (82) 1 , it is obvious that this system has at least two roots since Eq. (82) 1 can be rewritten as follows where F + 1 and F − 1 respectively denote the function F ± 1 associated with the signs + and − before the square root. We propose and discuss in AppendixE a strategy to choose the appropriate root.
In addition to this issue, we have to face to an unexpected problem related to a possible non positivity of the secant modulus µ (2) 0 in some cases. As an illustration, Fig. 2 depicts the variations 5 of this modulus as a function of time, in the case an elastically reinforced incompressible composite with an ideally elastoplastic matrix submitted to a cyclic isochoric macroscopic extension. It is observed that the value of secant modulus µ (2) 0 is sometimes negative, especially for t ∈ [14. 6; 22]. Of course, the occurrence of negative values for the secant modulus of a LCC does not comply with the classical framework underlying the homogenization theory of heterogeneous materials. Although a theoretical analysis on the possibility of negative shear modulus associated with a linear thermo-elastic fictive LCC is out of the scope of the present paper, we can however make the following remarks.
First, the LCC of energy w 0 for which µ (2) 0 is sometimes negative is a fictitious and not a real composite which is determined by mathematical stationarity conditions which do not take into account any constraints on the sign of its modulus. Second, in such a situation, Eq. (53) 2 does in principle no longer lead to the definition of an effective thermoelastic composite, whose properties and local field statistics would be evaluated for instance by the Hashin and Shtrikman equations provided in AppendixD. We nevertheless propose to formally make use of these relations in such cases and we observe that doing so, a solution of the nonlinear system Eq. (82) can be obtained. Finally and more importantly, it is observed that this solution is consistent not only in terms of macroscopic behavior but also local fields statistics with numerical predictions obtained with FFT or FE simulations or with other MFH approaches, see Section 4. This fact, especially the good agreement at the local scale of the field statistics, suggests that the use of linear MFH approaches -such as the Hashin and Shtrikman estimate employed in this work to homogenize the LCC of energy w 0 -when formally applied to isotropic thermoelastic LCCs with negative shear modulus is consistent in a sense that remains to be determined. Accordingly, the notion of LCC even defined with negative secant modulus -and therefore purely fictitious -can still remain a relevant tool in the framework of nonlinear homogenization.

Validation by comparaison to predictions of previous models
This section aims to compare the results obtained by the DIV approach to the ones obtained by either other MFH approaches of the literature or full-field FFT or FE simulations. For that, the DIV approach is applied to elasto-plastic incompressible composites made of isotropic elastic spherical inclusions, randomly and isotropically distributed in an elasto-plastic matrix. The latter can be either elastic ideally-plastic, or with linear kinematic or isotropic hardening. The applied loadings are cyclic, radial or non radial.

Loading histories
Two types of loading are considered. The first one has been considered by Lahellec and Suquet [16] and consists of an applied macroscopic strain E(t) composed of an isochoric extension E 33 (t) along the axial direction and of equal shears E 13 (t) along the two planes parallel to the axial direction In association with (88), two types of loading histories, presented in Fig. 3, are considered. The first one corresponds to a cyclic radial extension while the second one is fully non-radial since it involves a rotation of the macroscopic strain and therefore of the principal axes of the macroscopic stress because of the overall isotropy of the composite.
The second loading to be considered is a tensile loading It should be noted that the use of a tensile loading requires to adapt the algorithm presented in 3.4.1 since the latter is based on a macroscopic strain loading. For that, we proceed in the following way. Due to the isotropy of the composite, the macroscopic strain associated with the tensile stress (89) reads where α is an unknown scalar.
and shows that parameter α only depends on µ ε p ) defined by Eq. (82) 1 also depends on the macroscopic strain. Since the latter only depends on µ (2) 0 and η (2) ε p , finally F 1 still only depends on the two unknows µ (2) 0 and η (2) ε p of the nonlinear problem. A similar argument holds for the function F 2 defined by Eq. (82) 2 such that it becomes possible to derive the root of the nonlinear problem by the procedure fsolve and to compute the macroscopic tensile stress Σ 11 by means of Eq. (91) 1 when the macroscopic strain E 11 (t) is imposed.

Case of ideally-plastic matrix
In this section we consider the composite material having an elastic ideallyplastic matrix as studied by Lahellec and Suquet [16] and defined by the following material parameters except that now the inclusion and matrix phases in our simulation are incompressible In what follows, predictions of the DIV approach are compared with both the RVP formulation developed by Lahellec and Suquet [16] and Fast Fourier Transform (FFT) simulations carried out in the same work on a RVE comprised of 50 spherical inclusions randomly distributed in the matrix. Fig. 4a depicts the overall response of the composite under isochoric cyclic extension while the evolutions of the average stresses in the matrix and inclusion phases are shown in Figs. 4b and  4c, respectively. Lastly, Fig. 4d illustrates the evolution of the stress fluctuations in the matrix. On a whole, a very good agreement is observed between the predictions of the DIV and RVP formulations. When compared to full-field simulations, it is observed that the DIV approach, similarly to the RVP formulation, provides accurate predictions for the evolutions of both the macroscopic stress and the average stress within the matrix, even if the macroscopic response is slightly overestimated while the reverse is true for σ 33 (2) . Moreover, both DIV and RVP approaches capture the slight Bauschinger effect observed on the macroscopic response. However, although the DIV and RVP approaches qualitatively reproduce the trends of the FFT simulations, they strongly overestimate the average stress within the inclusion as well as the stress fluctuations within the matrix. For the evolution of σ 33 (1) , it is noted that the DIV approach during the plastic reloading fits better to the FFT simulations than the RVP approach. The reverse occurs if we consider the evolution of the stress fluctuations since the RVP formulation manages to reproduce the local minimum values of the FFT stress fluctuations at the beginning of the plastic reloadings (t ∼ 14, 5s and t ∼ 34, 5s) unlike the DIV approach which predicts almost zero fluctuations in this part of the curve. Lastly, as can be seen in Figs. 4a and 4b, the slopes associated with the elastic regime are slightly different for the DIV and RVP formulations. This is simply due to the difference in Pois-son's ratio used for the simulations, namely ν (1) = ν (2) = 0.4999 for DIV while ν (1) = ν (2) = 0.3636 for RVP.
At this stage, it is worth mentioning that, except for the elastic slope, the differences between both compressible (DIV) and incompressible (RVP and FFT) simulations are expected to be negligible or at least small since carried out on reinforced composites submitted to macroscopic loadings with zero (isochoric extension) or low stress (tensile loading) triaxiality as it is the case in our study. To ensure this point, FE simulations have been performed on two-phase periodic composites composed of spherical elastic particles embedded on a perfect cubic lattice in an elasto-plastic matrix, with or without hardening, submitted to a cyclic macroscopic isochoric extension and whose material properties are defined by Eq. (93), and by Eq. (94) or (96) in case of linear kinematic or isotropic hardening, respectively. Both compressible ν = 0.3636 and incompressible ν = 0.499999 cases have been considered. As expected, it is found that the difference between the compressible and incompressible evolutions of macroscopic and local responses during the plastic regime are negligible (less than 0.2%) for both the macroscopic axial stress and average axial stress within the matrix and small (less than 7%) for the average axial stress in the inclusion and the stress fluctuations within the matrix -which take into account all the components of the local deviatoric stress. Theses numerical observations justify the comparisons made on these macroscopic and local quantities between our DIV approach evaluated for incompressible composites and the RVP and FFT simulations evaluated for the same composites but with compressible phases.
As seen above, the DIV and RVP models are in a very close agreement on significant parts of the cyclic curves. A succinct theoretical comparison between both models can help to provide some explanations for this close agreement. Indeed, although both models rely on distinct LCCs, it is easily seen that they share some points in common such that the same definition of their respective secant viscosity (see Eq. (37)), the same evolution of the cumulative plastic strain p (see Eq. (38) 3 ) and the fact that both models satisfy the plastic yield criterion when applied to the second-order moment σ − X (r) (see Eq. (44)). These common features are mainly due to the fact that both DIV and RVP approaches make use of the variational procedure, two times for the DIV approach and one time for the RVP. Furthermore, although they rely on two different incremental variational principles, these principles are close since the one used by the RVP approach can in a first approximation be interpreted as a rate version of the one used by the DIV procedure and initially proposed by Lahellec and Suquet [13]. It should be noted that the comparisons reported in Fig. 4 have been carried out for an elastic contrast which is fairly small E (1) /E (2) = 2. In order to see the influence of this parameter on the macroscopic and local responses of the composite, additional simulations of the DIV model with elastic contrasts of 5, 10, 50 have been performed and compared with FE computations carried out on a periodic cubic cell made of a single spherical elastic inclusion embedded in an elastic ideally-plastic matrix. The numerical results reported in Figure 5 clearly show that the more the elastic contrast increases, the less the Bauschinger effect is captured by the DIV model. Furthermore, for high elastic contrast, it seems that the DIV approach fails to reproduce qualitatively the evolution of the stress fluctuations during the beginning of the plastic reloading.

39
For the non-radial loading associated with program 2, we report in Fig. 6 the evolution of the macroscopic axial and shear stresses w.r.t. the time. It is observed that the DIV approach is in good agreement with the RVP approach and FFT simulations even if it always slightly overestimates (resp. underestimates) the FFT data for the macroscopic axial stress (resp. for the macroscopic shear stress).

Matrix with linear kinematic hardening
We now consider an elasto-plastic matrix with a linear kinematic hardening characterized by the scalar a (2) introduced in Eq. (62). The composite parameters are the same as in Section 4.1.2 and the kinematic hardening parameter is set as in Lahellec and Suquet [16] a (2) = 300 MPa.
Predictions of the DIV approach are again compared with the FFT simulations carried out by Lahellec and Suquet [16] on a RVE comprised of 50 spherical inclusions randomly distributed in the matrix. For the case of a radial extension (program 1), Fig. 7 depicts the evolutions of the macroscopic axial strain, the average stresses over the phase as well as the stress fluctuations. The same trends as those already obtained in Fig. 4 for the reinforced composite with an ideally elasto-plastic matrix are observed except the fact that the axial stress average in the matrix now exhibits a linear kinematic work-hardening. All comments made in the previous section for the reinforced composites with an ideally plastic matrix still hold for the linear kinematic hardening case except that now the DIV approach provides more accurate predictions of both the average stress in the inclusion and the matrix fluctuations than for an ideally-plastic matrix. The variations of the back-stress, depicted in Fig. 8, show that its average in the matrix is in good agreement with the FFT simulations. However, its fluctuations in the matrix are strongly overestimated -a shift of about 100% is observed -even if they qualitatively follow the same trends as the FFT simulations.

Matrix with isotropic hardening
We now consider an elastoplastic matrix with an isotropic hardening characterized by the power law As a first illustration, we consider the reinforced composite dealt with by Lahellec and Suquet [16] whose material parameters are defined by Eq. (93) for the elastic properties of the phases and by the following parameters for the isotropic hardening β (2) = 100 MPa, γ (2) = 0.4 .
First, we apply the radial macroscopic cyclic extension associated with program 1 (see Fig. 3a). For such a loading, Fig. 9 depicts the evolutions of both the axial effective stress and the average stresses in each phase as well as the stress fluctuations in the matrix for the DIV approach and for both the RVP formulation and FFT simulations computed [16]. We observe the same trends as those obtained in Fig. 4 for an elastic ideally plastic matrix except that now the axial stress average in the matrix exhibits an isotropic hardening. All comments made for the reinforced composites with an elastic ideally-plastic matrix still hold for the isotropic hardening case. Fig. 10 depicts the evolutions of the macroscopic axial and shear stresses obtained for the non-radial loading associated with program 2 of Fig. 3b. It is observed that the DIV approach is in good agreement with the FFT simulations. However, the FFT data are always slightly overestimated for the macroscopic axial stress and underestimated for the macroscopic shear stress.
We now consider the Aluminium / Silicone Carbide Metal Matrix Composite (MMC) dealt with by Brassart et al. [6] for a cyclic tensile test and compare the

Discussion
It has been seen in the previous sections that both RVP and DIV formulations provide close and accurate predictions of the local and macroscopic responses. It would be interesting to try to determine the situations for which there is an advantage to use one or the other of these methods. The following remarks provide some insight about this issue. For compressible matrix, as reported in AppendixE, the DIV approach is not fully numerically stable since it encounters, at some scarce points, convergence problems due to a singularity of the function F 2 (see Eq. (82) 2 and Fig. E.18). This is not the case of the RVP approach which is numerically stable. For incompressible matrix, both the DIV and RVP approaches are numerically stable and both provide accurate results for elasto-plastic composites with or without local hardening. However, when the matrix exhibits kinematic hardening, it seems that the RVP approach is a bit less efficient numerically than the DIV procedure since the incorporation of the kinematic hardening within the RVP approach is made at the price of two supplementary unknowns H (r) 0 and X (r) n (see [16]) which increase the numerical cost of the method, unlike the DIV approach for which no supplementary unknown is required to handle the linear kinematic hardening.

Multiple cycles loadings
In this Section, we explore the predictions of the DIV model when several loading cycles are considered: program 1 on Fig. 3a is reproduced several times during the considered loading sequences. For that, we compute both the macroscopic and local responses of the reinforced composites considered in the previ- For an elastic ideally-plastic matrix or an elastoplastic matrix with linear kinematic hardening, it is found that the macroscopic and local responses are stabilized as soon as the first cyclic loading occurs. This is not the case for an elastoplastic matrix with isotropic hardening as it can be seen on Fig.  12 which reports for 10 cycles the evolution of the axial stress, the averages of the axial stress over the phase as well as the stress fluctuations in the matrix for both the DIV approach and FE periodic simulations. As in former sections, the FE simulations are carried out on a periodic cubic cell made of a single spherical elastic inclusion embedded in an elastic ideally-plastic matrix. On a whole, a close agreement is observed between the DIV approach and the FE simulations, especially for the macroscopic axial stress and for the average axial stress in the matrix. As observed on the macroscopic axial stress, the asymmetry characterising the Bauschinger effect increases continuously with the number of cycles, going from 24,64 Mpa for the first cycle to 40,6 Mpa for the tenth cycle. This evolution of the asymmetry characterising the Baushinger effect is accurately captured by the DIV approach. Lastly, it can be seen that the macroscopic and local responses tends to a limit cycle, thus showing that the DIV approach, in agreement with the FE simulations, predicts a plastic shakedown when the plastic matrix exhibits isotropic 46 hardening.

Combined isotropic and linear kinematic hardening
This section provides new data for composites made of spherical linear elastic particles randomly distributed in an elastoplastic matrix with combined linear kinematic and isotropic hardening. The DIV approach is compared to FE simulations that we carried out on one eight of a three dimensional periodic cubic cell made of a centered inclusion surrounded by the matrix. Conditions of symmetry and periodicity are applied on this cell. The FE simulations are performed with the software Cast3M which makes use of the following isotropic hardening law The material parameters of the composites are set by Eq. (93) for the elastic properties and by the following equation for the linear kinematic and isotropic hardening R max = 2.1 GPa, β = 0.26, a (2) = 100 MPa.
The reinforced composite is subjected to a cyclic macroscopic isochoric extension (see program 1 of Fig.3). The results are reported in Fig.13 which depicts the variations of the effective axial stress, the average axial stresses over the phase as well as the stress fluctuations in the matrix. It is observed that the DIV macroscopic response is in good agreement with the FE simulations and captures well the Baushinger effect during the second plastic loading (E 33 ≈ 0.085). However, the DIV model fails to accurately reproduce the slope characterizing the work-hardening zones and for this reason underestimates the effective stress with a maximum error of about 3%. The average stress in the matrix obtained for the DIV approach is in very close agreement with the FE simulations except at the beginning of the third plastic loading (E 33 ≈ −0.075) where it fails to reproduce the Baushinger effect. The variation of the average stress in the inclusion derived by the DIV formulation qualitatively exhibits the same trends as the one observed for the macroscopic response. However the shift between the slopes which characterize the work-hardening zones becomes much more significant than for the effective response. This result was predictable. Indeed, since the average stress in the matrix is accurately evaluated by the DIV approach, the slight shift observed for the macroscopic response between the slopes which characterize the work-hardening zones should be strongly amplified in the inclusion for low volume fractions. This shift is probably due to the anisotropy of the problem to be solved by FE since the geometry of the unit cell is cubic unlike the microstructure addressed by the DIV approach which is assumed isotropic by the use of the Hashin and Shrikman isotropic estimate. If we now turn to the evolution of the stress fluctuations, it is observed that they are in close agreement with the FE simulations unlike it was the case for the previous composites considered in this study -i.e. made of either an ideally-plastic matrix or an elasto-plastic matrix with linear kinematic or isotropic hardening -for which only the trends of the FFT simulations were qualitatively but not quantitatively reproduced by the DIV approach. Of course, this statement holds for the whole cyclic loading except during the transition from the elastic regime to the plastic one (t ≈ 14s and t ≈ 37s) for which the stress fluctuations obtained by the DIV approach are again close to zero.

Conclusion
This study is based on the incremental variational procedure initially proposed by Lahellec and Suquet [13] to predict the effective behavior of dissipative composites. Recently, in order to evaluate the effective behavior of elasto-viscoplastic composites without hardening as well as their field statistics, Agoras et al. [1] proposed an alternative formulation to the modified incremental variational principles introduced by Lahellec and Suquet [16]. In the present paper, we extend the work of Agoras et al. [1] to the context of an elasto-plastic matrix with combined isotropic and linear kinematic hardening. Furthermore, the proposed extension is established in primal form (based on w(ε) and ϕ(α)) while a dual formulation was considered in Agoras et al. [1]. To this end, similarly to the work of Agoras et al. [1] the variational procedure of Ponte Castañeda [21] is applied two times: first, to linearize the nonlinear phase behaviour, including hardening, and to approximate the local and effective behaviors of the nonlinear composite by those of a LCC with intraphase heterogeneous polarization. The second application of the variational procedure allows to determine the behavior of the LCC with intraphase heterogenous polarization by making use of another LCC with homogeneous polarization whose material coefficients are derived by the stationarity conditions induced by the second application of the variational procedure. This approach is applied to elastically reinforced two phases isotropic composites and results in a system of three nonlinear equations with three unknowns which is numerically solved by using the Levenberg Marquardt algorithm.
The DIV approach has been first applied to three different composites, made of either an ideally plastic matrix or an elastoplastic matrix with linear kinematic or isotropic hardening submitted to two different loading histories corresponding to a proportional cyclic macroscopic isochoric extension and a non-radial loading, respectively. The predictions of the DIV approach have been compared to the results obtained by the RVP approach and FFT simulations [16]. On a whole, for all the composites considered a good agreement between the DIV and RVP formulations and the FFT simulations is obtained both at the macroscopic and local scales. In particular, the Bauschinger effect is well reproduced when the elastic contrast between the phases is not too large. However, the DIV approach needs to be improved since the average stress in the inclusion is overestimated as well as the stress fluctuations in the matrix which in addition show some sig-nificant deviations from the numerically simulated data at the beginning of the plastic reloading. For that, it might be appropriate to use more accurate linearization schemes than the variational procedure such as the recent partially or fully stationary second-order estimates of Ponte Castañeda [25,26]. Then, as in Brassart et al. [6] simulations for a tensile cyclic loading applied on a elastically reinforced metal matrix composite have been carried out and show a good agreement between the DIV approach and both their FE simulations and incremental variational homogenization scheme.
Finally, new data on elastically reinforced composites either with combined isotropic and linear kinematic hardening or submitted to several loading cycles have been proposed. On a whole, a good agreement is observed between the DIV approach and the FE simulations we carried out, even if the evolution of the average stress in the inclusion needs probably to be improved. Since it is now possible to assess accurately the effective behavior as well as the field statistics of elasto-plastic composites with combined isotropic and linear kinematic hardening, a natural continuation of this study would consist to tackle the case of nonlinear kinematic hardening. Such a work is now in progress.
The quantity τ being a scalar, we can use the classical KKT optimality conditions to solve the optimization problem (A.2), such that Taking into account that domains Ω and Ω (r) are fixed -and therefore that the derivative of the integral is equal to the integral of the derivative -and making use of the chain-rule, Eq. (A.3) 1 rewrites as where ∂ denotes the variation w.r.t. α = (ε p , p). Since Eq. (A.4) holds for arbitraryα, one gets Making use of the variable change λ (r) = c (r) λ (r) and choosingα = α, we obtain the stationarity conditions (39).
AppendixB. Determination of expression (52) for the effective thermoelastic energyw L (E) Let us introduce the function V (r) defined as The stationarity is achieved forê 0 . By considering the above result, function V (r) can be rewritten as Note that this approximation is a lower (resp. upper) bound for L (r) 0 satisfying ∆L (r) ≥ 0 (resp. ∆L (r) ≤ 0). The effective behaviour is therefore given bỹ Let A be a fourth-order tensor. In this appendix, we seek to demonstrate the following identity where I 4 denotes the symetric fourth-order identity tensor. Moreover where I 8 is the eighth-order identity tensor defined as I 8 = I 4 ⊗ Noting that I 4 : B = B for any eighth-order tensor B, we have AppendixD. Effective behavior and field statistics for thermo-elastic linear composites Let us consider a N-phase thermoelastic composite whose per-phase freeenergy is given by According to Willis [30], its effective free-energy reads w 0 (ε) = 1 2 E :L 0 : E+τ 0 : 0 tensors denote the classical strain localization tensors. As shown for instance [23], the first and secondorder moments of the strain over the phase can be obtained by the following relationsε (r) = A where subscripts (1) and (2) denote the inclusion and matrix phases, respectively, and ∆L 0 = L 0 . In this study, the microstructural tensor P

AppendixE. Multiple solutions
To well understand the consequences induced by the non uniqueness issue raised in Section 3.4.3 (see Eq. (87)) and to define a strategy allowing to choose the appropriate root, we consider here the simplest situation which corresponds to the case of an elastic ideally-plastic matrix. For that, simulations are carried out on the same elastically reinforced composite considered in Section 4.1.2 and submitted to a cyclic isochoric macroscopic extension. As shown in Fig. E.14 which depicts the macroscopic response, a non physical jump is observed around E 33 ≈ 0. The same jump also appears for the first and second-order moments of the stress over the phase, and for the stress fluctuations in the matrix. This jump is due to the non uniqueness of the root of the equation F 1 µ (2) 0 , η (2) ε p = 0 with F 1 defined by (82) 1 . In fact, it is numerically observed that the function F 1 has several roots w.r.t variable µ (2) 0 at fixed η (2) ε p : two as expected for composites with incompressible phases (see Eq. (87)) but even four for compressible composites.
Case of incompressible composites. Let us first focus on composites with incompressible phases. For this kind of material, it has been numerically observed that the function F 1 has always two solutions as illustrated by Fig. E.15 which displays the variations of the function F 1 w.r.t. µ (2) 0 for fixed η (2) ε p at the value E 33 corresponding to the jump (see E.14). In order to get rid of the jump observed in Fig. E.14, we need to select the appropriate root. For that, we add in the algorithm the following constraints Σ 33 (t n+1 ) − Σ 33 (t n ) E 33 (t n+1 ) − E 33 (t n ) ≥ 0 and σ  which allow to enforce the continuity of both the macroscopic stress and the stress average over the matrix, and therefore to check if the root provided by the procedure fsolve of matlab is the appropriate one. If constraints (E.1) are satisfied, we can go the next time step. If not, we have to determine the other root of F 1 since the root initially provided by fsolve is not the good one. For that, we simply change the value of the starting point associated with the variable µ (2) 0 which is required by the procedure fsolve. As shown in Fig. E.16, this correction of the algorithm yields good results since the evolution of the macroscopic stress is now regular with no jump occurrence as before. Moreover, the additional numerical computing cost induced by this correction is weak: around 3s for a total computational time of 20s.
Case of compressible composites. We now consider the case of compressible composites. For that, we applied the DIV approach to the same composite considered in the previous paragraph -that is the elastically reinforced composite with an elastic ideally-plastic matrix defined by the set of material parameters (93) -except that now the Poisson's ratios of the phase are given by ν (1) = ν (2) = 0.3636. The loading is still a cyclic isochoric macroscopic extension. It is numerically observed that the solution of equation (2) 0 for fixed η (2) ε p at time t = 20s corresponding to E 33 = 0. The determination of the appropriate root becomes more difficult because of the occurence of four roots. To obtain the appropriate one, the same constraints (E.1) than the ones presented for incompressible composites are used. When these constraints are not satisfied by the solution initially provided by the procedure fsolve, we now have to determine the three other remaining roots and check which of them satisfies constraints (E.1). To determine the three remaining roots, we make use of a dichotomy method which, together with constraints (E.1), allow us in most of the cases to obtain the appropriate root. However, for some scarce points the algorithm fails to determine the appropriate root. This occurs because the function F 2 (see also Eq. (82) 2 ) presents for compressible composites a singularity, as illustrated on Fig. E.18 which reports the evolution of this function w.r.t. µ (2) 0 for fixed η (2) ε p at time t = 20s corresponding to E 33 = 0. The analytical expression of F 2 suggests that the singularity is probably generated by a pole. However, for the time being, we did