Spurious softening in the macroscopic response predicted by the additive tangent Mori–Tanaka scheme for elastic–viscoplastic composites

The Mori–Tanaka (MT) scheme is a well-established mean-field model that combines simplicity and good predictive capabilities. The additive tangent MT scheme is a popular variant of the method that is suitable for elastic–viscoplastic composites. This work is concerned with the analysis of some intrinsic features of the additive tangent MT scheme, in particular, of spurious softening in the macroscopic response that may be encountered when the Perzyna-type viscoplasticity model is used. The resulting non-monotonic macroscopic stress–strain response is clearly non-physical, but it also has a negative impact on the efficiency and robustness of the MT model when it is used as a local constitutive model in concurrent multiscale finite-element computations. As shown in the paper, the spurious softening is more pronounced when the so-called soft isotropization is employed to compute the viscoplastic Hill tensor, but it is also observed, although for a much narrower range of material parameters, in the case of the hard isotropization and when no isotropization is applied. Moreover, the softening is promoted at low strain rates, for high elastic contrast, and for high volume fractions of inclusions. Nevertheless, if the soft isotropization is avoided, the additive tangent MT scheme proves to be a feasible and computationally robust mean-field model that can be successfully employed in finite-element computations.


Introduction
When estimating effective properties of heterogeneous materials, mean-field models are known to offer fairly good predictive capabilities at low computational cost, as compared to more advanced methods, especially those based on the computational homogenization. Among the mean-field schemes applicable to two-phase composites of the matrix-inclusion geometry, the Mori-Tanaka (MT) method (Mori and Tanaka, 1973) is the most widely used. It makes this mean-field scheme the first-choice candidate for application in large-scale finite-element analyses.
The additive tangent MT model belongs to the third category and is potentially one of the best candidates for implementation as a material model at element Gauss points in finite-element calculations. The idea of the additive interaction law for heterogeneous elastic-viscoplastic materials was first proposed by Molinari et al. (1997) in the context of a self-consistent scheme for polycrystalline materials with isotropic overall properties. In Molinari (2002) the concept was adopted to a MT-type micro-macro transition scheme for composite materials. Basic advantages of the model are its appealing simplicity brought by the additive form of the interaction equation as well as the fact that, for an isotropic incompressible linear viscoelastic material, it delivers an exact solution to the inclusion problem (Mercier et al., 2005). It is worth mentioning here that  came up with a more involved interaction equation which delivers an exact solution of the inclusion problem also for compressible materials. In the context of non-linear elastic-viscoplastic two-phase composites, the additive tangent MT model has been recently successfully validated with respect to finite-element calculations in Czarnota et al. (2015) and Mercier et al. (2019) for uniaxial cyclic loadings. Probably one of the first finite-element implementations of the scheme with the Perzyna-type viscoplasticity model was performed by Msolli et al. (2016) using UMAT subroutine in ABAQUS/Standard. The tangent matrix is there computed using the small-perturbation technique. However, it seems that computational efficiency of the implementation and related issues have not been sufficiently discussed. It is worth noting that for the last decade the predictive capabilities of the additive interaction law have been also explored in conjunction with the self-consistent scheme and micromechanical modelling of polycrystalline materials (Wang et al., 2010;Abdul-Latif et al., 2018;Zecevic and Lebensohn, 2020;Jeong and Tomé, 2020;Girard et al., 2021).
The goal of the present study is to explore intrinsic features and possible simplifications of the additive tangent MT model that may affect computational efficiency of its finite-element implementation. In the analysis the Perzyna-type viscoplasticity model is assumed for the elastic-viscoplastic matrix phase, while inclusions are purely elastic.
Our recent works (Sadowski et al., 2017a,b) on the finite-element implementation of the incremental MT model of elastic-plastic composites have revealed convergence problems that turn out to result from discontinuities in the incremental finite-step response at the vicinity of the elastic-to-plastic transition and from the related abrupt change of the reference stiffness in the MT interaction equation. The problem has been overcome by proposing a sub-stepping strategy. Although similar discontinuities are not expected when using Perzyna-type model, which can be seen also as a viscoplastic regularization of rate-independent plasticity, possible issues related to the elastic-to-plastic transition, which may affect robustness of implementation, need to be identified. Such an issue, which is discussed in this work and, to our best knowledge, has been overlooked in the previous studies, is a spurious nonmonotonic macroscopic stress-strain response under strain-controlled proportional loading.
The original formulation employs the anisotropic viscoplastic tangent stiffness for which the respective Hill tensor involved in the additive interaction law needs to be calculated by performing numerical integration. Obviously, this has a negative impact on the model efficiency and robustness. Therefore, the isotropization techniques known for the rate-independent elastic-plastic tangent stiffness (Doghri and Ouaar, 2003) are of interest. It will be shown that application of some isotropization strategies may, however, lead to amplification of the spurious effects mentioned in the previous paragraph. Additionally, contrary to the elastic-plastic composites, isotropization seems not to improve the quality of predictions as compared to the anisotropic model (Czarnota et al., 2015). Therefore such a strategy must be selected with care.
The paper is organized as follows. The additive tangent MT model employing the Perzyna-type law for the matrix is presented in Section 2. Attention is paid to the anisotropic viscoplastic compliance and stiffness tensors as well as to the related Hill tensor. Possible strategies of model isotropization are also discussed. Section 3 illustrates selected effects resulting from the MT model, such as a non-monotonic (or softening) overall response of a composite material in the viscoplastic regime and a non-smooth elastic-to-plastic transition, not observed for a homogeneous elastic-viscoplastic material governed by the Perzyna model. Section 4 is devoted to a detailed analysis of the softening in the macroscopic response assuming the special case of proportional deviatoric loading for which the tensorial governing equations simplify to scalar ones, thus facilitating a more detailed analysis. In Section 5, predictions of the MT model with and without isotropization are compared to the results of unit-cell computations. Finally, in Section 6, an application of the additive tangent MT model in concurrent multiscale finite-element computations is presented in order to illustrate the impact of the spurious softening on the efficiency and robustness of the overall computational scheme.

Additive tangent Mori-Tanaka scheme
The MT model (Mori and Tanaka, 1973), a mean-field model originally developed for linear-elastic composites, relies on the Eshelby solution (Eshelby, 1957) to the problem of an ellipsoidal inclusion immersed in an infinite linearly elastic matrix. The Eshelby solution implies the following interaction equation (Hill, 1965), where the stress 1 and strain 1 in the inclusion are homogeneous, and those in the matrix, 0 and 0 , are the far-field quantities. The fourthorder tensor L e * is called the Hill tensor (or the interaction tensor) and depends on the elastic properties of the matrix and on the inclusion shape. The MT model for a linear-elastic composite material is then obtained by identifying 0 and 0 with the average stress and strain in the matrix, and the system of equations is closed by imposing the averaging scheme and the local constitutive equations of the phases.
Consider now a two-phase composite material with elastic spherical inclusions embedded in a Maxwell-type elastic-viscoplastic matrix. The local constitutive equations in the rate form thus reaḋ where the total strain in the matrix is additively decomposed into elastic and viscoplastic parts,̇0 =̇e 0 +̇v 0 , and L 0 and L 1 are the fourthorder elastic stiffness tensors. In general, the viscoplastic parṫv 0 of the strain rate in the matrix is a non-linear function of 0 , the current stress in the matrix. For the purpose of the present work, the viscoplastic part of matrix deformation is assumed incompressible, so the tensoṙv 0 is deviatoric and depends only on the matrix stress deviator. For the considered elasto-viscoplastic material behaviour, the following additive form of the interaction equation is used in the present work, Here, M e * is the inverse elastic Hill tensor, M e * = (L e * ) −1 , see Eq.
(1), which depends on the elastic properties of the matrix, M v * is the inverse viscoplastic Hill tensor, which depends on the viscoplastic compliance tensor of the matrix and will be specified later, and = − 1 3 ( tr ) are the stress deviators in the phases with denoting the second-order identity tensor. The additive law (3) was first formulated by Molinari (2002) and next the additive tangent MT model, adopted in the present study, was elaborated in full by Mercier and Molinari (2009). Note that the same interaction equation is obtained in the context of the MT averaging scheme when using the concept of sequential linearization due to Kowalczyk-Gajewska and Petryk (2011) (for details, see Section 2.2 and Variant II of the model in that paper).
The current overall strain̄and stress̄of the composite material are updated by averaging the corresponding local rate quantities,̄= where is the volume fraction of inclusions. The interaction equation (3), the averaging rule (4) and the constitutive Eqs.
(2) complemented by the constitutive law for the viscoplastic strain ratėv 0 form a complete set of equations that govern the overall response of the considered elastic-viscoplastic two-phase composite.

Perzyna-type viscoplastic model
As mentioned above, the matrix is assumed to be governed by an elastic-viscoplastic model, and the Perzyna-type overstress model (Perzyna, 1963) is here adopted for that purpose. For consistency with the subsequent derivations, the subscript '0' (referring to the matrix) is explicitly used in the constitutive equations below.
The elastic domain, 0 ≤ 0, is specified by the 2 (Huber-von Mises) yield function, wherêy 0 is the yield stress that may in general depend on the (isotropic) hardening parameter 0 , and ‖ 0 ‖ = √ 0 ⋅ 0 . Here and below, a centred dot denotes the scalar product. Evolution of v 0 is governed by the plastic flow rule, where the plastic multiplieṙ0 is an explicit function of the stress deviator 0 and is given bẏ Here, 0 is the relaxation time and is the rate-sensitivity parameter, 0 < ≤ 1. The hardening parameter 0 is defined as the accumulated plastic strain with the evolution laẇ0 =̇0.

Viscoplastic compliance and stiffness tensors
As a step towards introducing the inverse viscoplastic Hill tensor M v * in the interaction equation (3), let us define the viscoplastic compliance tensor M v 0 , wherėv 0 is specified by the flow rule (6). Applying the spectral decomposition (Rychlewski, 1995;Kowalczyk-Gajewska and Ostrowska-Maciejewska, 2009), the fourth-order tensor M v 0 , understood as a linear operator in the six-dimensional space of symmetric second-order tensors, can be expressed in the spectral form, in terms of the orthogonal projectors G and F, and the corresponding inverse Kelvin moduli (i.e., the eigenvalues in the spectral decomposition), wherė0 is specified by Eq. (7). In Eq. (10), I is the symmetrized fourthorder identity tensor, while I P and I D are the orthogonal projectors, respectively, onto the spherical and deviatoric spaces of second-order symmetric tensors, so that I = I P + G + F and we have GF = O, GG = G and FF = F. It follows that the inverse Kelvin modulus associated with the spherical part I P is equal to zero, as a result of plastic incompressibility. The viscoplastic stiffness tensor L v 0 is then defined as the inverse of the compliance tensor M v 0 . However, considering that M v 0 is positivesemidefinite, the inverse is here defined in the five-dimensional deviatoric space, such that The use of the orthogonal projectors G and F makes the calculations straightforward. Accordingly, L v 0 reads The Kelvin moduli 2 1 and 2 2 , where 1 < 2 for ≤ 1, are of multiplicity one and four, respectively.

Elastic and viscoplastic Hill tensors
The elastic Hill tensor L e * is defined in the standard manner (Hill, 1965) in terms of the elastic stiffness tensor of the matrix, L 0 , and of the polarization tensor P e , namely L e * = (P e ) −1 − L 0 , P e =P(L 0 ), whereP denotes, in a symbolic way, a function that delivers the polarization tensor for a given argument (here the stiffness tensor L 0 ) and for a specified inclusion shape (assumed spherical throughout this work), see e.g. Willis (1981). In the case of elastically isotropic matrix and spherical inclusions, the elastic Hill tensor L e * and its inverse M e * can be readily expressed in closed form (e.g., Hill, 1965;Kowalczyk-Gajewska, 2012).
Likewise, the viscoplastic Hill tensor is defined in terms of the viscoplastic stiffness tensor L v 0 , Eq. (12), Similarly to Eq. (12) the formula is restricted here to the fivedimensional deviatoric space, which is admissible due to the proposed form of the interaction equation. Recall that, for an incompressible material, the Hill tensor in the six-dimensional space is well defined, and it has a finite bulk modulus (Hutchinson, 1976). However, only its projection onto the deviatoric space, I D L v * I D , plays a role in the interaction equation (3). Accordingly, the viscoplastic compliance Hill tensor M v * is derived as an inverse of the Hill tensor L v * in the five-dimensional deviatoric space, The polarization tensor P v in Eq. (14) depends on the viscoplastic stiffness of the matrix and on the shape of inclusions. Generally, L v 0 is anisotropic even if the material itself is isotropic, and thus the derivation of P v =P(L v 0 ) is only possible via numerical integration (Willis, 1981). Note that, in the case of an incompressible solid, in most of available procedures a very small compressibility is introduced to L v 0 (e.g., Hutchinson, 1976;Bornert et al., 2001). Such an approach is also followed in the present work. An alternative solution was proposed in Lebensohn et al. (1998), where the tangent stiffness given by Eq. (12) was used to calculate the polarization tensor, however, an additional constraint was imposed on Green's function to account for viscoplastic incompressibility.
In order to reduce the computational cost, the subsequent simplification is commonly used so that P v,iso is determined in terms of the isotropized stiffness of the matrix, L v,iso 0 , and the polarization tensor is then isotropic as well. As in the case of the elastic-plastic tangent stiffness in rate-independent plasticity, isotropization of L v 0 may be performed in various ways (e.g., Doghri and Ouaar, 2003;Chaboche et al., 2005;Pierard and Doghri, 2006;Czarnota et al., 2015).
Referring to the specific form (12) of the viscoplastic stiffness tensor L v 0 , its isotropic counterpart can be defined as where is defined in terms of the moduli 1 and 2 , and is a weighting factor, 0 ≤ ≤ 1. There are two common choices of factor . Setting = 0, so that = 1 , corresponds to the so-called soft isotropization (recall that 1 < 2 ). On the other hand, the hard isotropization is obtained for = 4∕5 so that the resulting tensor L v,iso 0 is the isotropic part of L v 0 in the sense that Now, for an isotropic viscoplastic stiffness tensor L v,iso 0 , the polarization tensor P v,iso and its inverse, both restricted to the five-dimensional deviatoric space, are given in an explicit form, Finally, upon isotropization of L v 0 , the viscoplastic Hill tensor can be determined directly using Eq. (14). By expressing the Hill tensor entirely in terms of the isotropic stiffness L v,iso 0 , we have Note that, alternatively, the isotropic polarization tensor, Eq. (18), can be combined with the actual anisotropic viscoplastic stiffness L v 0 , thus yielding an anisotropic Hill tensor, However, in the case of the soft isotropization, the Hill tensor given by Eq. (20) is not necessarily positive definite as required by the theory (Walpole, 1981). On the other hand, it has been checked that for the hard isotropization the Hill tensor of the form (20) leads to an overly stiff response. Accordingly, the second option, Eq. (20), is not used in the detailed analysis reported below, and only the isotropic Hill tensor (19) is used whenever isotropization is employed. Note that, in the context of isotropization, the names 'soft' and 'hard' stem from the fact that there exists a monotonous relation between the Hill tensor L * and the corresponding stiffness tensor L in the sense that if L − L is positive definite then L * − L * is also positive definite. As a consequence, the stress-strain response calculated for the soft isotropization will be softer than the one predicted for the hard isotropization. Based on the same reasoning, for the hard isotropization, one may expect a softer response for the case (19) than for (20).

Illustration of selected effects resulting from the Mori-Tanaka model
The additive tangent MT model outlined in the previous section is now applied to illustrate some effects that are characteristic for this micromechanical scheme. Specifically, our aim here is to demonstrate that in some situations, depending on the material and loading parameters, a softening response is predicted by the model, which can be considered a spurious effect of the MT model when applied to elastic-viscoplastic materials. A detailed study of the softening and of the related effects is carried out in Section 4.
The material parameters used in the present illustrative simulations refer to a metal-matrix composite with spherical (ceramic) elastic inclusions. The Young's modulus and the Poisson's ratio of the matrix are adopted as 0 = 75 GPa and 0 = 0.3, respectively, and those of the inclusions as 1 = 400 GPa and 1 = 0.2. The yield stress of the matrix is assumed equal to y 0 = 75 MPa, and ideal plasticity (i.e. no strain hardening) is assumed so that̂y 0 ( 0 ) = y 0 , unless otherwise stated. The rate-sensitivity parameter is adopted as = 1, unless otherwise stated. Finally, the relaxation time 0 is left unspecified, since it is used to normalize the strain rate, so that the latter is always specified by the product 0̇̄. The computations have been carried out using a constant strain increment Δ̄= 10 −4 . It has been checked that the results are not visibly affected when the strain increment is further reduced.
In the simulations reported below, unless otherwise stated, the composite material is subjected to uniaxial tension, and the macroscopic response is reported in terms of the macroscopic tensile stress̄and strain̄. The prescribed constant macroscopic strain rate 0̇̄i s varied and so is the volume fraction of inclusions. The MT model combined with the soft isotropization, see Eq. (19) with = 1 , is used in most cases, since this model admits an analytical solution that is discussed in more detail in Section 4.4, however, other cases are also considered. Fig. 1 shows the tensile response of the (homogeneous) matrix material ( = 0) for = 1. In this reference case, the analytical solution can be easily obtained (e.g., de Souza Neto et al., 2008, Section 11.2.7). The response is linear elastic until the stress reaches the yield stress y 0 , and the corresponding strain is denoted by y 0 (both quantities are used to normalize the response in Fig. 1). The elastic-to-plastic transition is smooth (the slope of the stress-strain curve is continuous), and then the stress increases towards its asymptotic value ∞ 0 = (1 + 0̇0 ) y 0 , see Section 4.3. This asymptotic response is described by a simple exponential function.
Now, Fig. 2 shows the corresponding response predicted by the additive MT model for a composite material with the inclusion volume fraction = 0.1. The elastic response, including the macroscopic yield stress̄y and the corresponding macroscopic strain̄y, is governed by the classical MT model, cf. Eq. (1), and is thus rate-independent. For the adopted material parameters, the macroscopic yield stress̄y is higher than that of the matrix, while the macroscopic yield strain̄y is lower as a result of a higher effective elastic modulus.
The response in the plastic range is qualitatively different than that in the reference case of a homogeneous elastic-viscoplastic material (i.e. for = 0). As expected, the asymptotic macroscopic stress̄∞ increases with increasing strain rate 0̇̄. However, there is a range of strain rates for which̄∞ is lower than the macroscopic yield stress̄y, see the case of 0̇̄= 0.01 and 0.05 in Fig. 2. As a result, in some cases, the stressstrain response is non-monotonic. This is a spurious effect of the MT model, since in general such a response is not expected in reality (and is not predicted by more elaborate micromechanical schemes, see e.g. the finite-element unit-cell computations reported in Section 5). To the best of our knowledge, this feature of the additive tangent MT scheme has not been reported nor discussed to date, and hence its extended analysis is carried out in the present paper. Fig. 2 shows two types of the non-monotonic (or softening) response. The macroscopic stress may decrease monotonically in the plastic range, as in the case of 0̇̄= 0.01, or first increase and then decrease towards̄∞, as in the case of 0̇̄= 0.05 and 0.1. For higher strain rates, e.g. for 0̇̄= 0.5, the macroscopic stress-strain response is monotonic.
We also note that, unlike in the case of a homogeneous material ( = 0), the elastic-to-plastic transition is not smooth, which is clearly visible in the case of the lower strain rates. The kink on the stress-strain curve is not visible (but still present, see Section 4.4) for 0̇̄= 0.5. Fig. 3 illustrates the effect of the volume fraction on the macroscopic response for a fixed macroscopic strain rate 0̇̄= 0.05. Note that in panels (a) and (b) the response is normalized by the yield stress and strain of the matrix ( y 0 and y 0 , respectively), hence the effect of on the effective elastic modulus, on the macroscopic yield stress and on the asymptotic macroscopic stress is clearly visible (all those quantities increase with increasing ). In panel (c), the response is normalized by the macroscopic yield stress and strain (̄y and̄y, respectively), hence the elastic response appears independent of and the effect of on the asymptotic macroscopic stress in reversed. It can be seen in Fig. 3 that the softening is more pronounced for higher volume fractions . We also note that for = 0.05 the asymptotic macroscopic stress is higher P. Sadowski et al.    than the macroscopic yield stress,̄∞ >̄y, still the response exhibits softening after the macroscopic stress reaches a maximum at̄∕̄y ≈ 1.1, see the horizontal dashed line in Fig. 3c.
To check whether the softening response is not associated with the rate-sensitivity parameter set to = 1 in the preceding simulations, the effect of parameter is examined in Fig. 4 for = 0.1 and for three selected strain rates. It can be seen that the macroscopic response is affected by the rate-sensitivity parameter and by the strain rate 0̇ī n a complex manner, but the general observation is that softening may occur also for < 1.
The effect of strain hardening of the matrix material is illustrated in Fig. 5. Specifically, a linear isotropic hardening is considered such that the current yield stress, cf. Eq. (5), is defined aŝy 0 ( 0 ) = y 0 + 0 0 , where 0 is the hardening modulus (recall that ideal plasticity with 0 = 0 has been assumed so far). As expected, depending on the hardening modulus 0 and on the magnitude of the softening in the respective case with 0 = 0, the hardening effects in the matrix may fully compensate the softening in the macroscopic response so that a monotonic response is obtained, as in the case of 0 = 0.05 0 for 0̇̄= 0.1 in Fig. 5b, or only partially reduce the magnitude of softening, as in the remaining cases depicted in Fig. 5a,b. In general, the softening during the initial transient period (if present) is followed by a hardening overall response thus a non-monotonic up-down-up response is obtained for 0 > 0.
Finally, we examine the impact of the isotropization method applied to the viscoplastic stiffness tensor L v 0 that is used to compute the Hill tensor L v * . Here, the case of the 'soft' and 'hard' isotropization is considered, the second one in two variants corresponding to Eqs. (19) and (20), as well as the case of no isotropization (denoted as 'anisotropic'). In the latter case, the predictions vary depending on the assumed deformation process, and here two limit cases of uniaxial tension and pure shear are considered. The results presented in Fig. 6 may suggest that softening occurs only for the soft isotropization and in the case of pure shear with no isotropization. However, a detailed analysis carried out in the next section shows that there exists a (small) range of material parameters for which softening may be observed also in other cases. This is illustrated in Fig. 7 in which the material parameters have been adopted such that softening, even if relatively weak, is observed in all cases, except for the hard isotropization combined with Eq. (20). The results reported in Figs. 6 and 7 confirm also that the hard isotropization combined with Eq. (20) yields an overly stiff response, and hence this case is not discussed in the subsequent sections.
In the illustrative cases reported so far, only the macroscopic response has been examined. To provide more insight into the effects under consideration, Fig. 8 presents the evolution of the local stresses in the phases as a function of the macroscopic strain for the case studied in Fig. 6. From Fig. 8 and also from other observations that are not reported here, it follows that the softening in the macroscopic response is accompanied by a decrease of the equivalent stress in the (elastic) inclusions during the initial stage of deformation soon after the elasticto-plastic transition. The related effects are more pronounced for the soft isotropization, but are also observed in other cases. This seems to be related to the release of the internal stresses that accumulate in the inclusions during elastic loading as a result of the elastic mismatch between the phases. Note that the equivalent stress in the matrix is practically not affected by isotropization or by the Lode angle in the anisotropic case, cf. Fig. 8a. Very small differences, illustrated in the insets in Fig. 8a, occur only in the vicinity of the elastic-to-plastic transition.
The main conclusion of the preliminary analysis carried out above is that the additive tangent MT model for elastic-viscoplastic composites may lead to a softening macroscopic response. For the rate-sensitivity parameter = 1, this effect is more pronounced at lower strain rates and at higher volume fractions of inclusions. However, the results shown in Fig. 4 indicate that for = 0.01 the effect is observed also at higher strain rates. We also note that the associated stress drop is often very abrupt and thus would not be compensated by the strain hardening in the matrix. It has also been observed that the macroscopic response predicted by the additive tangent MT model may exhibit a non-smooth elastic-to-plastic transition. Recall that the elastic-to-plastic transition is smooth in the case of a homogeneous elastic-viscoplastic material.

Analysis of the softening in the macroscopic response
In this section, we carry out a detailed analysis of the softening response resulting from the additive tangent MT model, which has been illustrated in the previous section. To this end, the special case of proportional deviatoric loading is considered, which results in scalar governing equations, Section 4.1. The elastic-to-plastic transition is discussed in Section 4.2. In Section 4.3, general expressions for the asymptotic macroscopic stress̄∞ are derived and the conditions under which̄∞ is lower than the macroscopic yield stress̄y are analysed. Finally, in Section 4.4, an analytical solution is derived for the case of soft isotropization and = 1, which facilitates a complete characterization of the softening response in this case.

Proportional deviatoric loading
Consider now the case of proportional deviatoric loading such that the composite material is subjected, in a strain-controlled process, to a deviatoric macroscopic strain, P. Sadowski et al.  where the tensor is a constant unit deviator, tr = 0, ‖ ‖ = 1, and the scalar factor̄=̄( ) is a prescribed monotonic function of time (̇̄( ) > 0), thus specifying the macroscopic strain̄. On account of plastic incompressibility and elastic and plastic isotropy of the constituent phases, it can be shown that, upon isotropization of the viscoplastic Hill tensor, the microscopic and macroscopic response is also deviatoric and proportional, see also Sadowski et al. (2017b), so that in particular, and likewise for the rates. Here,̄, and are scalar coefficients specifying the magnitudes of the respective deviatoric tensors. By construction,̄and correspond to the respective equivalent uniaxial stresses.
As a result, the tensorial governing equations of the MT model, Section 2, simplify to the respective scalar forms. The constitutive Eqs. (2) and (6) take the following form, where are the elastic shear moduli such that 2 = ⋅ L anḋ0 is governed by Eq. (7). The interaction equation (3) readṡ where e * and v * are scalar coefficients of the form 2 e * = ⋅ L e * = 0 (8 0 + 9 0 ) 3(2 0 + 0 ) , Here, 0 stands for the bulk modulus of the matrix, and is defined in Eq. (17). Coefficient e * corresponds to an elastically isotropic matrix and Eq. (25) 1 is a standard result (e.g., Hill, 1965). Coefficient v * of the form (25) 2 corresponds to the isotropic viscoplastic Hill tensor (17). Note that a closed-form expression is not available in the case of an anisotropic viscoplastic Hill tensor (14) which must be then evaluated numerically (Ponte Castañeda, 1996).
Finally, the averaging rule (4) is applied to the scalar quantities,̄= The set of Eqs. (23)-(26) can now be used to obtain the overall response of the considered elastic-viscoplastic two-phase composite under proportional deviatoric loading,̄= wherēis the overall elastic shear modulus, whilė̄v denotes the macroscopic plastic strain rate. Note thaṫ̄v is the overall inelastic strain rate, thuṡ̄v =̇̄−̇̄e, wherė̄e =̇̄∕(3̄), and P. Sadowski et al.  is not equal to the simple average of the local quantities in the form (1 − )̇v 0 .
Remark 1. Let us stress that, when the anisotropic viscoplastic Hill tensor is used in the interaction equation (3), the local and overall stresses and strains remain coaxial (they share common principal directions) for the considered process (21), but are not necessarily proportional, as exposed by Eqs. (22). It can be proved that Eqs. (22) hold also for the anisotropic viscoplastic Hill tensor when the isochoric tension/compression or pure shear processes are assumed and do not hold for the combination of these two processes. Strictly speaking, Eqs. (22) where denotes the Lode angle. Since by Eq. (15) the polarization tensor and thus the Hill tensor share their symmetry group with the stiffness tensor L v 0 , it is seen that Eq. (29) holds when ( ) represents isochoric tension/compression ( = 0 or ∕3) or pure shear ( = ∕6). In these two cases, L v * is transversely isotropic or of tetragonal symmetry, respectively, so ( ) is an eigenstate and 2 v * ( ) is the corresponding Kelvin modulus of any fourth-order tensor of the respective symmetry group. Note that the value of v * is in general different in these two cases. In other instances, L v * is orthotropic (with the orthotropy axes coaxial with the principal directions of ) and , even when being the eigenstate of L v 0 , is not necessarily an eigenstate of the polarization tensor P v and of the Hill tensor L v * . A direct consequence of the anisotropy of the Hill tensor is the influence of the third invariant on the composite equivalent strain-equivalent stress response, despite the local constitutive laws (2) and (5)-(7) do not involve such dependency. Note that the equality of eigensubspaces of anisotropic P v and L v 0 was erroneously claimed by Ponte Castañeda (1996) and Ponte Castañeda and Suquet (1997). This error was next corrected in (Nebozhyn and Ponte Castañeda, 1999) where it was shown that the assumption of such an equality provides only an approximation of P v .

Elastic-to-plastic transition
In the following, a monotonic loading process is considered witḣ̄> 0 and with the initial state specified bȳ= 0 and v 0 = 0. Initially the response is elastic, and the corresponding macroscopic response,̄= 3̄̄, cf. Eq. (27), results from the classical MT model, cf. Eq. (1). The macroscopic stress̄y at the instant of the elastic-to-plastic transition is then obtained by enforcing the yield condition 0 = 0 and reads P. Sadowski et al. As illustrated in Section 3, the macroscopic response may be nonsmooth at the instant of the elastic-to-plastic transition. From Eq. (27), it follows that the condition for smoothness is that the plastic strain rate vanishes at this instant, i.e.,̇̄v | | 0 =0 = 0. To examine this condition further, let us consider the viscoplastic response in the matrix at the elastic-to-plastic transition, i.e., the first branch in Eq. (7) at 0 = 0 (so when locally in the matriẋ0 = 0). The corresponding macroscopic plastic strain ratė̄v = √ 3∕2̇̄v can be found in a closed form for the isotropized viscoplastic Hill tensor and readṡ̄v It can be checked thaṫ̄v | | 0 =0 ≥ 0. Furthermore, we havė̄v | | 0 =0 > 0 only for the soft isotropization with = 1 (provided > 0 and 1 ≠ 0 , to exclude these two trivial cases). Indeed, for the soft isotropization ( = 0), by evaluating v * in Eq. (25) 2 using Eqs. (7), (11) and (17), we have which for = 1 simplifies to 2 v,sof t * = 0̂y 0 .
In the case of the hard isotropization ( = 4∕5), we have Accordingly, we have v * → ∞ for 0 → 0 + in all cases except for the soft isotropization with = 1, when v * is finite, cf. Eq. (33). Specification (31) still holds for the anisotropic Hill tensor under the condition (29), i.e., for a uniaxial tension/compression or pure shear process with v * = v * ( ). For other values of , at the instant of the elastic-to-plastic transition the proportionality (22) is lost, sincė̄v ≠ √ 3∕2̇̄v . As mentioned above, iḟ̄v | | 0 =0 = 0 then the macroscopic response is smooth at the elastic-to-plastic transition, cf. Eq. (27). The nonsmooth response, as observed in some cases in the sample simulations in Section 3, is thus limited to the case of the soft isotropization with = 1. Otherwise, the macroscopic response is smooth and thus shares this feature with the response of a homogeneous elastic-viscoplastic material.

Asymptotic macroscopic stress̄∞
In this section, we consider the case of ideal plasticity (no hardening,̂y 0 = y 0 ) and we examine the asymptotic macroscopic stress ∞ under a constant macroscopic strain rate, which, as shown below, can be determined in a closed form. We then provide the conditions in which̄∞ is lower than the macroscopic yield stress̄y, which is a sufficient (but not a necessary) condition for the softening in the macroscopic stress-strain response.
First, we consider the case when the Hill tensor is anisotropic, however, the condition (29) is fulfilled so that the proportionality relations (22) remain valid. Under a constant macroscopic strain rate 0̇̄, an asymptotic steady state is reached for̄→ ∞. In the steady state, the stresses are constant, thuṡ0 =̇1 =̇̄= 0, so that we havė0 −̇v 0 = 0 anḋ1 = 0, cf. Eq. (23). Accordingly, the interaction equation (24) and the averaging rule (26) simplify tȯ where v * ( ) is defined by Eq. (25) 2 with = ( ), and is the Lode angle equal to 0 (or ∕3) and ∕6 for uniaxial tension (or compression) and pure shear processes, respectively. In such case, where For other values of , the condition (29) is not fulfilled and in such case Eq. (36) provides only the work-conjugate equivalent stress such that̄∞ ⋅̇̄=̄∞( )̇̄, however,̄∞ ≠ √ 2∕3̄∞( ) . Recall that a closedform expression for v * ( ) is not available so that it must be evaluated numerically.
In agreement with Eqs. (36)-(37), a higher value of̄∞ for uniaxial tension than for pure shear was found by Nebozhyn and Ponte Castañeda (1999) using second-order variational estimates for twophase viscoplastic composites, and also noted by Ponte Castañeda and Suquet (2001) in FFT-based full-field simulations reported therein. Note that, although the second stress moments were used by Nebozhyn and Ponte Castañeda (1999) to obtain the estimates of the flow stress, the discussed sensitivity to the Lode angle was a direct consequence of the anisotropy of the tangent stiffness. This feature is shared by the present model employing only the first moments. The final predictions will thus be qualitatively the same for both approaches, while the use of the second-order moments may lead only to quantitative differences in the predicted values.
Next, we consider the case of the isotropized viscoplastic Hill tensor. In this case, proportionality conditions (22) hold, and the value of v * defined by Eq. (25) 2 is independent of . By combining Eq. (35) with Eqs. (7) and (26) 2 , and with the formula for v * , Eq. (32) or (34) depending on the isotropization variant, the asymptotic macroscopic stress is obtained in the following form ∞,sof t = y 0 for the soft isotropization, and in the following form for the hard isotropization. In the limit case of 0̇̄→ 0, the asymptotic stresses simplify tō∞ ,sof t = y 0 and̄∞ ,hard = y 0 (1 + 6 5 ). It is noted that, upon isotropization, the response predicted by the MT model, including the asymptotic stress̄∞, does not depend on the Lode angle , while it does depend on in the anisotropic case. This can be considered a limitation of the isotropized additive tangent MT model. However, in the anisotropic case, the predicted impact of on ∞ is in most cases mild. The maximum difference between the limit cases of = 0 (tension) and = ∕6 (shear) increases with decreasing strain rate and with decreasing . The extreme case, independently of , is obtained for the strain rate approaching zero, when the ratio of ∞ for shear and for tension is equal to 1∕(1 + 0.923 ), according to Eq. (36) evaluated numerically. For a realistic range of volume fractions and for the normalized strain rate 0̇̄g reater than 0.01, the maximum difference does not exceed 10%, and typically it is visibly smaller.
The asymptotic macroscopic stress̄∞ can now be compared to the macroscopic yield stress̄y, Eq. (30), and the condition̄∞ <̄y implies the following condition for the macroscopic strain rate, where the coefficients and depend on the material properties and on the isotropization variant, viz.
An interpretation of the condition in Eq. (40) is provided in Fig. 10 which shows three domains in the space parameterized by 1 ∕ 0 and 0̇̄. The shaded domain in the lower-right corner indicates the combination of parameters for which̄∞ <̄y, as specified by the condition in Eq. (40). The diagram in Fig. 10 corresponds to the case of the soft isotropization with = 1 in which the line separating the other two domains can be found in terms of the analytical solution discussed in Section 4.4. Fig. 11 shows the lines of 0̇̄= ⟨ ⟩ 1∕ corresponding to the soft isotropization for selected values of and . It can be seen that the domain of̄∞ <̄y increases with increasing , which is consistent with the results reported in Fig. 3. On the other hand, the domain of̄∞ <̄y significantly decreases with decreasing (note the logarithmic scale in Fig. 11b). The corresponding diagrams for the hard isotropization (not shown for brevity) indicate that the domain of̄∞ <̄y is then significantly smaller than in the case of the soft isotropization. This is consistent with the observation that the condition > 0 is satisfied only for relatively high ratios 1 ∕ 0 and only for a limited range of 0 , specifically for 0 > 13∕35 ≈ 0.37 (which includes the case of elastically incompressible matrix, 0 = 0.5), see also Fig. 9.
In Fig. 12, the lines bounding the respective domains̄∞ <̄y are shown for the case of no isotropization, cf. Eqs. (30) and (36), for = 0 (tension) and = ∕3 (pure shear). As a reference, the respective lines corresponding to the soft isotropization, cf. Fig. 11a, are also included in Fig. 12. It follows that the range of parameters for which the condition̄∞ <̄y is met is smaller in the anisotropic case than in the case of the soft isotropization and it is smaller for tension (or compression) than for pure shear, which are two limit cases. This is consistent with the sample results reported in Figs. 6 and 7, where the softening is weaker (or even absent) in the anisotropic case, as compared to the soft isotropization. Concluding, we recall that the condition̄∞ <̄y is a sufficient, but not a necessary condition for the softening in the macroscopic response. As illustrated in Figs. 2-5, the response may exhibit a maximum followed by a softening branch towards an asymptotic valuē∞ that is higher than the macroscopic yield stress̄y (this case is studied in detail below for the soft isotropization with = 1). Nevertheless, the analysis carried out above and the condition in Eq. (40) in particular give an indication which range of material parameters may lead to softening.

Analytical solution for the soft isotropization with = 1
In the case of the soft isotropization with = 1 and no hardening (̂y 0 = y 0 ), the viscoplastic Hill tensor L v,iso * is constant and so is the scalar coefficient v * = v,sof t * = 0 y 0 ∕2, cf. Eq. (33). As a result, the governing equations are then linear and can be reduced to a system of linear ordinary differential equations with constant coefficients, viz.
where the components of matrix and vector are given in Appendix A. Note that matrix depends on the elastic properties of the phases, on the yield stress y 0 and on the volume fraction , while vector additionally depends on the strain rate 0̇̄, but does not depend on the yield stress y 0 . The analytical solution can now be easily found by enforcing the initial conditions at the instant of the elastic-to-plastic transition (at = 0), namely 0 (0) = y 0 and̄(0) =̄y. In particular, the macroscopic stress̄is found in the following form, where I and II are two relaxation times (equal to the eigenvalues of ) and I and II are the corresponding integration constants (the respective formulae are quite involved and are not provided here). The response exhibits thus a qualitative difference with respect to the case of a homogeneous elasto-viscoplastic material in which only one relaxation time is involved, namely 0 = ∞ 0 + I − ∕ I . It can be checked that, in the limit of = 0, we have I = 0 y 0 ∕(3 0 ) and II = 0 so that the analytical solution (45) reduces to the analytical solution for a homogeneous material (e.g., de Souza Neto et al., 2008, Section 11.2.7).
The analytical solution can now be used to examine the situation when the macroscopic stress increases, reaches a maximum, and then decreases towards the asymptotic stress̄∞. Such a situation occurs iḟ̄> 0 for = 0 anḋ̄< 0 for → ∞. If the relaxation times are ordered such that II > I > 0, this holds for II > 0 and II ∕ II < − I ∕ I . The maximum, if exists, occurs at = max , where max is determined from the conditioṅ̄= 0 thus yielding max = I II II − I log (  − I II   II I ) .
The limit case of no softening corresponds to II → 0 and the corresponding macroscopic strain rate 0̇̄c an be determined numerically by solving the non-linear equation II = 0 for given , 1 ∕ 0 , 0 and y 0 . Sample results are provided in Fig. 13. For the interpretation, refer to Fig. 10, where the line corresponding to the condition II = 0 separates the domain of no softening (upper-left corner) from the domain of max >̄∞ >̄y. Comparing Fig. 13 to Fig. 11a, it can be seen that P. Sadowski et al.   the range of parameters in which softening occurs is significantly larger than the domain of̄∞ <̄y, as also illustrated in Fig. 10.
Referring again to Fig. 10 we notice that the zones in which softening occurs (the shaded zones in Fig. 10) appear only for 1 ∕ 0 > 1, i.e. only when the inclusions are elastically stiffer than the matrix. Indeed, it can be checked that the two conditions sof t = 0 and II = 0 at 0̇̄= 0 are satisfied exactly for 1 ∕ 0 = 1 (recall that those two conditions are relevant for the soft isotropization, the latter additionally for = 1). In the case of the hard isotropization, the condition hard = 0 can only be satisfied for 1 ∕ 0 > 1, see Fig. 9. The respective zones of̄∞ <̄y do not exist for 0 < 13∕35 ≈ 0.37, see Fig. 9, and otherwise are relatively small.

Spherical unit cell
In this section, the predictions of the MT model in the various versions discussed above are compared to the results obtained for a spherical unit cell illustrated in Fig. 14. Motivated by the composite sphere model (Hashin, 1962), see also Kursa et al. (2018), the spherical unit cell is composed of a single elastic spherical inclusion located in the centre of the unit cell and surrounded by an elastic-viscoplastic matrix. The unit cell is subjected to isochoric tension along the -direction so that the problem is axisymmetric, and the corresponding finite-element computations are carried out for one quarter of the cross-section, as sketched in Fig. 14, with adequate symmetry conditions enforced along = 0 and = 0. Two canonical types of boundary conditions are applied. In the case of the linear displacement boundary conditions (LDBC), the displacement on the boundary is prescribed as =̄, wherēis the macroscopic strain, cf. Fig. 14b. In the case of the uniform traction boundary conditions (UTBC), the traction on the boundary, =̄, is derived from the macroscopic stress̄, where denotes the unit outer normal, cf. Fig. 14c. However, since the material response is here ratedependent, a special treatment is needed in the UTBC case in order to ensure that the macroscopic strain ratė̄is constant and prescribed to a desired value. Specifically, the UTBC have been prescribed by enforcing the so-called minimal kinematic boundary conditions (Mesarovic and Padbidri, 2005) which are actually equivalent to the UTBC, except that the macroscopic strain̄is prescribed and the macroscopic stressī s obtained as a result of the micro-macro transition scheme (Wojciechowski and Lefik, 2016). Accordingly, the same macroscopic strain rate (assumed constant) can be applied in both cases (LDBC and UTBC) thus facilitating a meaningful comparison of the results. As demonstrated by Ostoja-Starzewski (2006), for a proportional loading process for which the prescribed macroscopic strain̄is given by Eq. (21), the LDBC and UTBC lead to, respectively, an upper and lower bound on the effective properties. It should be noted that, strictly speaking, P. Sadowski et al.  the notion of bounds applies here to the particular configuration of the unit cell adopted and to the related composite microstructure in which a specific gradation of composite sphere sizes down to an infinitesimal size is assumed so that the whole representative volume is filled (Hashin, 1962).
The unit-cell problem is solved using the finite element method, and the standard incremental formulation of the small-strain Perzynatype elasto-viscoplasticity is employed (de Souza Neto et al., 2008). The standard details of the finite-element implementation are omitted here for brevity.
Material parameters are adopted equal to those used in Section 3. In particular, as before, no strain hardening is considered so that the transient response during the elastic-to-plastic transition is apparent. Figs. 15 and 16 show the results obtained for = 0.1 and for three representative macroscopic strain rates 0̇̄= 0.01, 0.1 and 1. The ratesensitivity parameter is adopted equal to = 1 in Fig. 15 and to = 0.1 in Fig. 16. In each case, the MT results are reported for the soft and hard isotropization as well as without isotropization with the viscoplastic modulus v * ( ) corresponding to uniaxial tension and pure shear. The case of shear is provided for completeness only, considering that the unit cell in finite-element calculations is subjected to tension. Note that, when the soft or hard isotropization is applied in the MT model, the resulting response is the same independently of the Lode angle . As mentioned above, the results corresponding to the finiteelement unit-cell model provide a kind of upper and lower bound when the LDBC and UTBC, respectively, are used, hence a shading has been added in Figs. 15 and 16 to indicate the possible intermediate values.
Black markers superimposed on the stress-strain curves indicate the limits of the elastic response. In the case of the unit-cell models, the markers correspond to the first instant at which the state changes from elastic to plastic in at least one Gauss point.
The results reported in Figs. 15 and 16 again illustrate the spurious softening response resulting from the MT model with the soft isotropization and with no isotropization (pure shear case). Importantly, for = 0.1, the softening is observed for all strain rates. Of course, the softening is not observed in the case of the unit-cell model. Moreover, the plastic zone in the matrix grows gradually, which leads to a more pronounced apparent hardening in the unit-cell models. Clearly, this effect is missing in the MT model in which uniform deformation is assumed in the matrix. In all cases, in view of no strain hardening, the macroscopic stress tends to its asymptotic value that depends on the macroscopic strain rate.
Interestingly, for = 0.1, the results corresponding to the three macroscopic strain rates shown in Fig. 16 are qualitatively very similar, only the stress level increases with increasing strain rate. This observation can be explained as follows. Qualitative differences between the stress-strain curves obtained for the anisotropic and isotropized viscoplastic Hill tensor result from the degree of anisotropy of the tensor L v 0 , Eq. (12), which can be quantified by the factor = ( 2 − 1 )∕ 2 , where 0 ≤ < 1. For the Perzyna viscoplasticity model, this factor can be specified in terms of and the equivalent strain rate in the matrix as follows: P. Sadowski et al. It is easily verified that, for a decreasing , the strain rate has a secondary impact on thus in such cases the shape of stress-strain curves is governed mainly by . (Czarnota et al., 2015) In this section, comparison is made with respect to the results of the finite-element computations carried out by Czarnota et al. (2015) for a periodic unit cell with 30 spherical inclusions. Here, both the matrix and the inclusions are elastic-viscoplastic and are governed by the Perzyna model with linear isotropic hardening. The material parameters can be found in Table 3 in Czarnota et al. (2015), and a constant overall strain rate is prescribed, 0̇̄= 3.33. The periodic unit cell is shown in Fig. 17a.

Periodic unit cell with 30 inclusions
In the description of the MT model in Section 2, the inclusions have been assumed elastic. Extension of the model to the case of elastic-viscoplastic inclusions, employed in the present example, is immediate and amounts to replacing the elastic constitutive equation of the inclusions, Eq. (2) 2 , by the elastic-viscoplastic one. Eqs. (3) and (4) and the definitions of M e * and M v * remain unchanged. Anyway, the effect of plastic deformations in the inclusions is here relatively small because their yield stress is significantly higher than that of the matrix, see Czarnota et al. (2015).
Predictions of the MT model with the soft and hard isotropization and with no isotropization are compared to the finite-element results of Czarnota et al. (2015) in Fig. 17. Two volume fractions of inclusions have been considered, = 0.1 and 0.25, and in both cases the soft isotropization leads to the softening in the stress-strain response, which is observed at each instant of the elastic-to-plastic transition. During the first loading, the softening is partially mitigated by the strain hardening in the matrix. However, during the subsequent unloading and reloading, the softening is visibly stronger, which is due to the higher internal stresses that build up in the microstructure. The softening is not observed in the other two cases. These two cases show also a very good agreement with the finite-element results, while the agreement is significantly worse in the case of the soft isotropization. Note that the very good agreement between the predictions of the additive tangent MT model employing the anisotropic viscoplastic Hill tensor and the finite-element results has already been demonstrated by Czarnota et al. (2015). The only difference in the present treatment of the additive tangent MT model with respect to that of Czarnota et al. (2015) is the time integration scheme, namely the implicit backward-Euler scheme has been used here, while the explicit forward-Euler scheme was used by Czarnota et al. (2015).

Mori-Tanaka model in multiscale finite-element computations
In this section, concurrent multiscale finite-element computations are carried out in order to illustrate the impact of the spurious softening on the efficiency and robustness of the overall computational scheme. To this end, the additive tangent MT model has been implemented in a finite-element code as a material model at the Gauss-point level. The respective incremental formulation has been obtained by applying the implicit backward-Euler integration scheme to the rate equations presented in Section 2. More details are provided in Appendix B, see also Sadowski et al. (2017a).
Sample computations are here carried out for a rectangular plate of the length = 20 mm and cross-section of 10 × 1 mm with a hole of the diameter of 5 mm, see Fig. 18a. The plate is loaded in tension in the displacement-control mode with a constant elongation ratė∕ . A fine mesh of about 370,000 eight-node hexahedral elements is used, which leads to nearly 1.2 million degrees of freedom (note that, for a better visibility, Fig. 18a shows a coarse mesh with the element size increased 4 times with respect to the actual mesh used in the computations). Accordingly, with eight Gauss points per element, the finite-element model involves about 3 million Gauss points at which the incremental equations of the MT model are solved at each global Newton iteration. It is thus a severe test of the robustness of the model both in terms of the local convergence at each Gauss point and in terms of the global performance. Recall that the spurious effects of the MT model are more pronounced at very small strain rates, and thus, for a problem with a highly non-uniform deformation, the chance of encountering such effects increases with an increasing size of the problem.
The composite material, governed by the MT model, is composed of elastic spherical inclusions and elastic-viscoplastic matrix. The elastic properties are specified by the Young's modulus and Poisson's ratio of the matrix, 0 = 75 GPa and 0 = 0.3, and inclusions, 1 = 400 GPa and 1 = 0.2. Perzyna-type viscoplastic model is adopted for the matrix with the isotropic hardening specified bŷ y 0 ( 0 ) = y 0 + 0 0 + ∞ 0 (1 − − 0 0 ), with the following hardening parameters:   is halved if convergence is not achieved in 12 iterations. The number of time steps needed to complete the simulation can thus be used as an indicator of the robustness of the model. Figs. 19 and 20 show the load-elongation curves obtained for the soft and hard isotropization, respectively. As a reference, the results obtained for a homogeneous material are also included in the figures. These results are labelled as = 0. The respective computations have been carried out using the MT model with = 0 and also using the classical model of elasto-viscoplasticity. The obtained responses are identical and so is the time-stepping history, which illustrates the robustness of the present implementation of the MT model.
In the case of the soft isotropization, the computations could not be completed for = 1 due to convergence problems caused by the spurious softening discussed in the preceding sections, see the case of = 0.2 in Fig. 19a. For = 0.1, the effect of the softening is weaker and only leads to an increase of the number of time steps needed to complete the simulation. This is more pronounced in the case of the low loading rate 0̇∕ = 0.01 which required 18 time steps, as compared to 5-11 time steps needed in all the remaining cases. Recall that the soft isotropization with = 1 leads not only to the softening in the macroscopic response (observed for = 0.1 as well), but also to a nonsmooth elastic-to-plastic transition, see Section 4.2. The latter effect may adversely affect the convergence behaviour and may explain the significant difference in the robustness of the model in the two cases.
As discussed in the preceding sections, the spurious softening is practically not observed in the case of the hard isotropization. Accordingly, no difficulties are encountered during the simulations, and the performance of the MT model is then comparable to that of the simple elastic-viscoplastic Perzyna-type model.

Conclusion
We have explored intrinsic features and possible simplifications of the additive tangent Mori-Tanaka model that affect computational efficiency of its finite-element implementation as a material model at element Gauss points. In the analysis, the Perzyna-type viscoplasticity model has been assumed for the elastic-viscoplastic matrix phase, while the inclusions are purely elastic. In particular, an issue related to the elastic-to-(visco)plastic transition has been identified, which significantly affects robustness of the model. Specifically, a spurious nonmonotonic macroscopic stress-strain response under strain-controlled proportional loading is predicted for some sets of material parameters and loading conditions. This type of behaviour is not observed in the reference case of a homogeneous elastic-viscoplastic material nor in more elaborate micromechanical schemes, e.g., in finite-element unit-cell computations. Such non-monotonic macroscopic response predicted by the MT model is often associated with an abrupt stress drop which cannot be compensated by strain hardening of the constituent  phases. Moreover, a non-smooth elastic-to-plastic transition has been observed at the macro level, while a smooth transition is expected, as it is characteristic for elastic-viscoplastic materials.
Since the MT model is a mean-field model, the stress and strain fields are assumed to be homogeneous within the constituent phases, and stress concentrations at inclusions are thus disregarded. This feature of mean-field models seems to be responsible for the spurious effects under consideration. In particular, in the MT model the elasticto-plastic transition occurs at a single well-defined instant when the yield criterion is met in the (homogeneous) matrix. This is in contrast to the gradual growth of plastic zones, which is captured by the fullfield models, such as those employed in Section 5. Accordingly, the macroscopic yield stress may be significantly overestimated by the MT model, see e.g. Figs. 15 and 16. When combined with the release of the internal stresses that accumulate in the elastic regime due to the elastic mismatch between the phases, this may lead to the softening in the macroscopic response. Actually, related effects may possibly be characteristic also for other mean-field models. Note that a nonmonotonic response similar to that observed in this work has been recently reported by Zecevic and Lebensohn (2020), see their Fig. 3, in the context of a self-consistent model which apparently also employs the additive tangent interaction law of Molinari et al. (1997), Molinari (2002). The mentioned deficiencies of the mean-field approach may possibly be mitigated, at least partially, by resorting to mean-field models employing second-order moments (e.g., Ponte Castañeda, 1996;Lahellec and Suquet, 2007), but these models are significantly more involved, which makes their application in concurrent finite-element computations problematic.
The original formulation of the additive tangent MT model employs an anisotropic viscoplastic Hill tensor derived from an anisotropic viscoplastic compliance tensor. Calculation of the Hill tensor is thus necessarily performed by means of numerical integration, which negatively affects efficiency and robustness of the computational scheme. Accordingly, the possible strategies of model isotropization have been discussed. As it has been demonstrated, such strategies must be employed with care, since some of them may lead to amplification of the spurious effects mentioned above. Specifically, such negative impact has been found to be particularly pronounced in the case of the soft isotropization. The related effects are more severe for the ratesensitivity parameter = 1 and for relatively low strain rates. On the contrary, for the hard isotropization, this effect is only observed in a narrow range of material parameters. Additionally, the latter isotropization method has been shown to provide a significantly better agreement with the finite-element unit-cell computations reported by Czarnota et al. (2015). It should be noted that, contrary to the MT model for rate-independent elastic-plastic composites, in the case of elastic-viscoplastic composites, isotropization seems not to improve the quality of predictions, as compared to the model with no isotropization.
Considering that the soft isotropization amplifies the spurious softening effects and delivers less accurate predictions, a general conclusion is that it should be avoided in the context of MT-based modelling of elastic-viscoplastic composites. On the other hand, the spurious effects are practically eliminated by the hard isotropization. At the same time, the hard isotropization has been found to provide satisfactory predictions in a realistic case studied in Section 5.2. The stiffening introduced by the hard isotropization is more visible in the range of parameters studied in Section 5.1, however, this example considers a less realistic case of no strain hardening at relatively low strain rates at which the stiffening effect is more pronounced. Importantly, the corresponding MT model is significantly cheaper and more robust than that with no isotropization. Hence the additive tangent MT model with the hard isotropization can be considered a reasonable compromise between accuracy on one hand and efficiency and robustness on the other hand.
The additive tangent MT model has been finally employed in concurrent multiscale finite-element computations. The respective incremental formulation has been obtained by applying the fully implicit backward-Euler integration scheme to the rate equations of the model. The results obtained for a sample boundary value problem illustrate a significant impact of the spurious softening on the efficiency and robustness of the overall computational scheme. On the other hand, in the softening-free cases, the present finite-element implementation of the MT model performs very well with the robustness comparable to the simple elastic-viscoplastic model. This has been tested for a relatively large-scale problem with the number of degrees of freedom exceeding one million.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The present finite-element implementation of the additive tangent MT model relies on the implicit backward-Euler integration scheme applied to the rate equations presented in Section 2. In the timediscrete formulation, the governing equations are enforced at time instant = +1 , while the rates are approximated by the respective finite differences, e.g.,̇= Δ ∕Δ , where Δ = +1 − and Δ = +1 − . The quantities pertaining to the previous time instant = (those with the superscript ) are known from the solution of the problem at the previous time step. The resulting non-linear algebraic equations are then solved using an iterative-subiterative scheme.
At the lowest level of constitutive equations of the individual phases, the incremental equations of Perzyna-type elasto-viscoplasticity are solved using the standard predictor-corrector algorithm (de Souza Neto et al., 2008, Section 11.6.4). The resulting stress-update algorithm is consistently linearized so that it can be effectively incorporated in the upper-level iterative scheme in which the interaction equation (3) is solved.
Application of the backward-Euler scheme to the interaction equation (3) leads to its following time-discrete form, where the viscoplastic compliance Hill tensor M v, +1 * is evaluated in terms of the current stress +1 0 . By exploiting the averaging rule (4) 1 , the local strains +1 can be expressed in terms of the overall strain +1 , considered to be known, is an unknown tensor. As a result, the timediscrete interaction equation (B.1) has been formulated as a non-linear equation with an unknown +1 , while the current stresses +1 depend on +1 through the stress-update algorithm described above with the local strains given by Eq. (B.2).
The Newton method is then applied to solve Eq. (B.2) and the macroscopic stress is obtained by averaging the local stresses,̄+ 1 = (1 − ) +1 0 + +1 1 . The corresponding macroscopic stress-update algorithm is consistently linearized, thus yielding the macroscopic tangent stiffness operator. This tangent operator is needed at the structural level at which the equilibrium equations are solved, again using the Newton method. The overall computational scheme is thus a doubly-nested iterative-subiterative Newton scheme, while consistent linearization at each level leads to a computationally efficient code exhibiting a quadratic convergence rate.
The computational treatment is here similar to that developed by Sadowski et al. (2017a) for the case of the rate-independent elastoplasticity. The overall structure is essentially the same, the differences are in the interaction equation (B.1), here resulting from the additive tangent MT scheme, and in the constitutive equations of the phases, here corresponding to Perzyna-type viscoplasticity. The reader is referred to Sadowski et al. (2017a) for the details of the computer implementation that is based on the automatic differentiation (AD) technique (Korelc and Wriggers, 2016).