We propose a micromechanics-based mean-field homogenization sch eme for the viscoelastic- viscoplastic particulate-reinforced composite which is applicab le to predict its mechanical response under complex loading conditions. We apply a formulati on based on an incrementally affine scheme by using algorithmic tangent operators, while ada ptively adjusting the strain of each constituent at every step of the loading process to ensure the consistency of the accumulated strain state and the concentration tensor. We name the method adaptive incrementally affine method. Despite mathematically rigorous de rivation, the method has some errors in plastic deformation regime. We propose an assumption for better prediction which dropping out the affine strain and affine stress in adaptive sc heme. We show that the adaptive incrementally affine scheme is able to predict the viscoelastic response very well. Still, it is inevitable that the plastic deformation of the composite is ini tiated earlier than our mean-field theoretical prediction because o f the local stress concentratio n near the particulate. Hence, we propose a yield reduction method that enforces the earlier init iation of the plastic deformation in the matrix phase when obtaining an effective mechanical resp onse. We show that the predictions from the adaptive incrementally affine scheme adjus ted with the yield reduction match well with various numerical simulations on particulate-re inforced composites considering viscoelastic, elasti c-viscoplastic, and viscoelasti c-viscoplastic matrices under uniaxial, cyclic, and bi-axial loadings. 3 1. Introduction Rate-dependent viscoelastic-viscoplastic behavior occurs in man y polymers and polymer composites. The viscoelastic-viscoplastic materials sho w not only a rate-dependent stress-strain response but also t he various characteristic resp onses such as stress relaxation under constant strain, creep with constant stress, and residual strain under unloading to zero stress. This paper aims to propose a micromechanics-based mean- field homogenization scheme so called adaptive incrementally affine method that can predict the effective mechanical behavior of composite consisting of the viscoelastic-viscoplast ic materials. For the linear elastic response, various mean-field homogenizat ion models have been studied, including the Eshelby me thod (Eshelby, 1957), Mori-Tan aka method (Mori and Tanaka, 1973), self-consistent method (Hershey, 1954; Hill, 1965), diff erential method (McLaughlin, 1977; Norris, 1985), and double inclusion method (Hori and Nema t-Nasser, 1993). The homogenization models for the linear elasticity are first exten ded to thermo-elasticity (Camacho et al., 1990; Lielens, 1999; Pierard et al., 2004), an d then to nonlinear responses involving viscosity or plasticity based on the similarity betwe en the thermo-elastic constitutive equation and the linearized equa tions concerning viscoelastic, plastic, and viscoplastic responses. The homogenization methods for the viscoelastic comp osites were developed based on the correspondence principle (Friebel et al., 2006; Hashin, 1965, 1970; Laws and McLaughlin, 1978). In this line of studies, the local constitut ive law in the time domain is converted to the Laplace-Carson domain before the application o f the homogenization scheme. Then, the effective behavior of composite in the time domain is calculated from the inverse Laplace transform. This framework is extended to study the elas to-viscoplastic composites (Masson and Zaoui, 1999; Pierard, 2006; Pierard and Doghri, 200 6a; Pierard et al., 2007). Unfortunately, the inversion of Laplace transform is mathematic ally complex even for the very 4 simple elasto-viscoplastic model, and its numerical implementat ion is complex and tricky for practical application. Hence, most recent studies have attempted to apply the homogeni zation scheme directly in the time domain wit hout relying on the Laplace-Cars on transform. An incremental method with a self-consistent scheme is proposed (Hill, 1965), Mori-Tanaka approach is applied with secant modulus (Tandon and Weng, 1988), and a simi lar approach with tangent modulus is suggested (Doghri and Ouaar, 2003). However, the pre diction from the aforementioned incremental methods turned out to be significant ly stiffer than the effective mechanical response of the composites measured in experiments o r computed from elaborate finite element analyses (FEA). To resolve the issue, an increme ntal secant method relying on fictitious secant moduli is suggested (El Ghezal et al., 2019; Wu et al., 2013a, b). However, the proposed secant method can only be applied to monotonic loading or proportional loading. Additionally, several other methods have been proposed includin g incremental variational principle (Lahellec and Suquet, 2007a, b, c), interaction law ( Mercier and Molinari, 2009; Molinari, 2002; Molinari et al., 1997), asymptotic homogenizati on (Frank and Brockman, 2001), and cell model (Kim and Muliana, 2010). For the sake of accuracy, some researches consider second-order moments reflecting the variation of field (Doghri et al., 2011; Wu et al., 2017; Wu et al., 2015). However, none of the aforementioned wor ks offered the solution for resolving both the stiff predic tion and limited loading issues. Recently, the incrementally affine method is devised as an alte rnative for resolving both problems (Doghri et al., 201 0; Miled et al., 2013). The pr imary advantage of the incrementally affine method is that it is applicable under mult i-axial loading as well. However, the predictions from the increm entally affine method also turne d out to be significantly stiffer than the experiments and the FEA results. To soften its predict ion, a heuristic isotropization 5 method was suggested that extrac ts the isotropic projection of anisotropic tangent operator for the homogenization procedure (Castañeda, 1996; Pierard and Dogh ri, 2006b). The regularization method which reduces the dependency of time incr ement has also been shown to soften the prediction as we ll (Doghri et al., 2010; Miled et al., 2011). However, both isotropization and regularizati on are empirical methods without rigorous mathematical reasoning or physical interpret ation (Chaboche et al., 2005). I t is found that, even with the isotropization and regularization methods, predictions from the incrementally affine method are still stiff. In the present study, we discussed the origin of the stiff pred iction from the incrementally affine method and suggested the adaptive incremen tally affine homogenization scheme that can be applied to ana lyze the viscoelastic-viscopla stic composites under complex loading conditions. We revealed that the inconsistency of accum ulated stress caused the incorrect prediction of effective stiffness of the composites, and proposed a method to adaptively adjust the strain of each constituent at every step of the loading process to ensure the consistency of the accumulate d strain state and the concent ration tensor (Kim, 2020). We found that the adaptive incrementally affine method is able to correctly predict the viscoelastic response, while it cannot correc tly describe the viscoplastic r esponse despite mathematically rigorous derivation. So, we proposed an assumption for better p rediction which dropping out the affine strain and affine stress in adaptive scheme. Althoug h non-uniform stress distribution exists within the matrix near the particle-matrix interface (he nce, the plastic deformation initiated locally near the particle), the s train/stress is dete rmined from the average strain/stress over the entire matrix region in the mean-field homogenization. So, we proposed a yield reduction method to tune the prediction on the viscoplastic res ponse, and validated the prediction from the homogenization scheme against the mechanica l response of viscoelastic- 6 viscoplastic composites under various loading conditions. The a daptive incrementally affine method has the following advantages. - Derivation without isotr opization or regularization - Applicability to the cyclic a nd multi-axial loading condition s - Accurate prediction of the m echanical response of composites with viscoelastic, elasto-plastic, elasto-viscoplastic, and viscoela stic-viscoplastic properties. This paper is organized as follows. Constitutive equations for modeling viscoelastic- viscoplastic behavior are summa rized in Section 2.1 Prony serie s is used for the linear viscoelastic model, and hardening function and viscoplastic fun ction are introduced for the viscoplastic model. The return mapping algorithm is presented i n Section 2.2. The adaptive incrementally affine method is described in Section 3. Rigorous derivation of adaptive incrementally affine method is presented in Section 3.1. Adapti ve incrementally affine method with an assumption for better prediction is introduced in Secti on 3.2. Sprit of yield reduction method is introduced in Section 3.3. Prediction of the adaptive incrementally affine method is validated against the finite element method (FEM) results, and the results are discussed in Section 4. Conclusions are drawn in Section 5. Boldface letters indicate second- or fourth-order tensor and no n-bold face letters indicate scalar values. (:) ref ers to double contraction and ( ⊗) refers to the dyadic product. 𝒂: 𝒃 ൌ 𝑎 ௜௝𝑏௝௜; ሺ𝑨:𝒃ሻ௜௝ൌ𝐴௜௝௞௟𝑏௟௞ሺ𝒂⊗𝒃 ሻ௜௝௞௟ൌ𝑎 ௜௝𝑏௞௟ 2. Material modeling and numerical method 2.1. Viscoelastic-viscoplastic modeling 7 We introduce the linear viscoela stic-viscoplastic model adopted in the present study, following the description of previous studies (Miled et al., 20 13; Miled et al., 2011). In the viscoelastic-viscoplastic const itutive model, the total strain appears as the sum of the viscoelastic strain ( 𝜺𝒗𝒆) and viscoplastic strain ( 𝜺𝒗𝒑) as follows, 𝜺ൌ𝜺𝒗𝒆൅𝜺𝒗𝒑. (1) In linear viscoelastic modeling, Cauchy stress 𝝈ሺ𝑡ሻ is obtained by using the Boltzmann’s hereditary integral (Boltzma nn, 1878) which include strain hist ory over [ െ∞, 𝑡] as 𝝈ሺ𝑡ሻൌන 𝑪 ሺ𝑡െ𝜏 ሻ:𝜕𝜺𝒗𝒆 𝜕𝜏𝜕𝜏௧ ିஶ. (2) Here, fourth-order stiffness tensor 𝑪ሺ𝑡ሻ for isotropic mater ials is defined as 𝑪ሺ𝑡ሻൌ2 𝐺 ሺ𝑡ሻ𝑰ௗ௘௩൅3 𝐾 ሺ𝑡ሻ𝑰௩௢௟ (3) 𝑰௩௢௟ൌ1 3𝟏⊗𝟏 ; 𝑰ௗ௘௩ൌ𝑰െ𝑰௩௢௟; where 𝟏 is second-order identity tensor and 𝑰 fourth-order identity tensor. 𝑰௩௢௟ and 𝑰ௗ௘௩ are volumetric and deviatoric p art of the fourth-order identity tensor, respectively. 𝐺ሺ𝑡ሻ and 𝐾ሺ𝑡ሻ are elastic shear modulus and bulk modulus, respectively. The shear modulus and bulk modulus can be expressed in the f orm of prony series as follows , 𝐺ሺ𝑡ሻൌ𝐺 ஶ൅෍𝐺 ௜exp ൬െ𝑡 𝑔௜൰ூ ௜ୀଵ (4) 𝐾ሺ𝑡ሻൌ𝐾 ஶ൅෍𝐾 ௝exp ቆെ𝑡 𝑘௝ቇ௃ ௝ୀଵ (5) where 𝐺ஶ a n d 𝐾ஶ are long-term shear modulus and bulk modulus, 𝐺௜ and 𝐾௝ a r e relaxation moduli, and 𝑔௜ and 𝑘௝ are corresponding relaxation times. Substituting the Eq. (3) into Eq. (2), stress can be decomposed into a deviatoric part, 𝒔ሺ𝑡ሻ, and volumetric part, 𝜎ுሺ𝑡ሻ as follows, 8 𝒔ሺ𝑡ሻൌ2 𝐺 ஶ𝝃𝒗𝒆ሺ𝑡ሻ൅෍𝒔 𝒊ሺ𝑡ሻூ ௜ୀଵ (6) 𝜎ுሺ𝑡ሻൌ3 𝐾 ஶ𝜀ு௩௘ሺ𝑡ሻ൅෍𝜎 ு௝ሺ𝑡ሻ௃ ௝ୀଵ (7) where 𝝈 ሺ𝑡ሻൌ𝒔ሺ𝑡ሻ൅𝜎 ுሺ𝑡ሻ𝟏; 𝜺𝒗𝒆ሺ𝑡ሻൌ𝝃𝒗𝒆ሺ𝑡ሻ൅𝜀ு௩௘ሺ𝑡ሻ𝟏 where 𝝃𝒗𝒆ሺ𝑡ሻ and 𝜀ு௩௘ሺ𝑡ሻ are the deviatoric part and volum etric part of viscoelastic st rain, respectively. Here, viscous stress 𝒔𝒊ሺ𝑡ሻ and 𝜎ு௝ሺ𝑡ሻ are defined as follows, 𝒔𝒊ሺ𝑡ሻ≡2 𝐺 ௜exp ൬െ𝑡 𝑔௜൰නe x p ൬𝜏 𝑔௜൰𝜕𝝃𝒗𝒆 𝜕𝜏d𝜏௧ ିஶ (8) 𝜎ு௝ሺ𝑡ሻ ≡ 3𝐾 ௝exp ቆെ𝑡 𝑘௝ቇන e x pቆ𝜏 𝑘௝ቇ𝜕𝜀ு௩௘ 𝜕𝜏d𝜏௧ ିஶ (9) One can rewrite Eqs. (8) and (9) to represent the update of vis cous stress values within a single time step ൫𝑡 ∈ ሺ𝑡௡,𝑡௡ାଵሻ൯ as follows, 𝒔𝒊ሺ𝑡௡ାଵሻൌe x p ൬െ∆𝑡 𝑔௜൰𝒔𝒊ሺ𝑡௡ሻ൅2 𝐺 ௜නe x p ൬𝜏െ𝑡 ௡ାଵ 𝑔௜൰𝜕𝝃𝒗𝒆 𝜕𝜏d𝜏௧೙శభ ௧೙ (10) 𝜎ு௝ሺ𝑡௡ାଵሻൌe x p ቆ െ∆𝑡 𝑘௝ቇ𝜎 ு௝ሺ𝑡௡ሻ൅3 𝐾 ௝නe x p ቆ𝜏െ𝑡 ௡ାଵ 𝑘௝ቇ𝜕𝜀ு௩௘ 𝜕𝜏d𝜏௧೙శభ ௧೙ (11) For the numerical integration of Eq. (10, 11), one can either a pproximate the strain rates ቀడ𝝃𝒗𝒆 డఛ,డఌಹೡ೐ డఛቁ within the time step ሺ𝑡௡,𝑡௡ାଵሻ to be constant (Frank and Brockman, 2001), or apply the mid-point rule (Simo and Hughes, 2006). We use the fi rst scheme for this study because the first method is rather insensitive to the time incr ement, ∆𝑡, under constant viscoelastic strain rate ( 𝜀ሶ௩௘) condition (Miled et al., 2011). Then, the viscous stresses in Eq. (10, 11) can be arranged as follows, 𝒔𝒊ሺ𝑡௡ାଵሻൌe x p ൬െ∆𝑡 𝑔௜൰𝒔𝒊ሺ𝑡௡ሻ൅2 𝐺 ෨௜∆𝝃𝒗𝒆 (12) 𝜎ு௝ሺ𝑡௡ାଵሻൌe x p ቆ െ∆𝑡 𝑘௝ቇ𝜎 ு௝ሺ𝑡௡ሻ൅3 𝐾 ෩௝∆𝜀ு௩௘ (13) where 𝐺 ෨௜ሺ∆𝑡ሻ≡𝐺 ௜൤1 െ exp ൬െ∆𝑡 𝑔௜൰൨𝑔௜ ∆𝑡, (14) 9 𝐾෩௝ሺ∆𝑡ሻ≡𝐾 ௝ቈ1െe x pቆ െ∆𝑡 𝑘௝ቇ቉𝑘௝ ∆𝑡. (15) Given the viscous stress, str ain, and strain increment at 𝑡௡, deviatoric stress and volumetric stress at 𝑡௡ାଵ are expressed as 𝒔ሺ𝑡௡ାଵሻൌ2 𝐺 ஶ𝝃𝒗𝒆ሺ𝑡௡ሻ൅෍e x p ൬െ∆𝑡 𝑔௜൰𝒔𝒊ሺ𝑡௡ሻூ ௜ୀଵ൅2 𝐺 ෨∆𝝃𝒗𝒆 (16) 𝜎ுሺ𝑡௡ାଵሻൌ3 𝐾 ஶ𝜀ு௩௘ሺ𝑡௡ሻ൅෍e x pቆ െ∆𝑡 𝑘௝ቇ𝜎 ு௝ሺ𝑡௡ሻ௃ ௝ୀଵ൅3 𝐾 ෩∆𝜀ு௩௘. (17) Here, incremental relaxation moduli 𝐺෨ and 𝐾෩ are defined as follows, 𝐺෨ሺ∆𝑡ሻ≡𝐺 ஶ൅෍𝐺 ෨௜ூ ௜ୀଵ (18) 𝐾෩ሺ∆𝑡ሻ≡𝐾 ஶ൅෍𝐾 ෩௝௃ ௝ୀଵ. (19) Having introduced the viscoelasti c part, we turn our attention to the viscoplastic response. Here, we adopted the classical 𝐽ଶ viscoplastic model with isotropic hardening, rate- dependent viscoplastic function. Yield function 𝑓 is defined as follows, 𝑓൫𝜎௘௤,𝑝 ,𝜺 ሶ൯≡𝜎 ௘௤െቀ 𝜎 ௬ሺ𝜺ሶሻ൅𝑅 ሺ𝑝ሻቁw h e r e 𝜎 ௘௤ൌඨ3 2𝒔: 𝒔 (20) where 𝜎௘௤ is von mises equivalent stress, 𝑝 is accumulated plastic strain, 𝜺ሶ is time derivative of strain, 𝜎௬ is yield stress, and 𝑅ሺ𝑝ሻ is hardening stress. The accumulated plastic strain ሺ𝑝ሻ is obtained by integrating the accumulated plastic strain rate ሺ𝑝ሶሻ over time as follows, 𝑝≡න𝑝 ሶሺ𝜏ሻd𝜏௧ ଴where 𝑝 ሶൌඨ2 3𝜺ሶ𝒗𝒑:𝜺ሶ𝒗𝒑 (21) For hardening response beyond the yield, isotropic hardening la w is adopted with the power- law as, 10 𝑅ሺ𝑝ሻൌ𝑘 𝑝௡ (22) where 𝑘 is hardening modulus and 𝑛 is hardening exponent. We note that the formulation can be readily extended to model nonlinear isotropic hardening and kinematic hardening (Doghri et al., 2010). Pla stic strain rate ሺ𝜺ሶ𝒗𝒑ሻ follows the plastic flow rule as follows, 𝜺ሶ𝒗𝒑ൌ𝑝ሶ𝜕𝑓 𝜕𝝈ൌ𝑝ሶ𝑵 (23) 𝑵≡𝜕𝑓 𝜕𝝈ൌ3 2𝒔 𝜎௘௤ (24) Here, the accumulated plastic strain rate ሺ𝑝ሶሻ is determined by a viscoplastic function, 𝑔௩൫𝜎௘௤,𝑝 ,𝜺 ሶ൯, which is given in terms of von mises equivalent stress, accum ulated plastic strain, and strain rate as follows, 𝑝ሶൌቊ𝑔௩൫𝜎௘௤,𝑝 ,𝜺 ሶ൯i f 𝑓൐0 0i f 𝑓൑0 . (25) We note that the methodology prop osed in this study can be appl ied to various forms of viscoplastic function ( 𝑔௩) such as Norton’s power law. In the numerical examples in the present study, we adopted the widely use d viscoplastic function with a power-law, 𝑔௩൫𝜎௘௤,𝑝 ൯ ൌ ቐ𝜅ቆ𝑓 𝜎௬൅ 𝑅ሺ𝑝ሻቇ௠ if𝑓൐0 0i f 𝑓൑0 (26) where 𝜅 and 𝑚 are viscoplastic modulus a nd exponent, respectively. 2.2. Return mapping algorithm for vi scoelastic-viscoplastic mat erials For the numerical tracking of the inelastic response, we applie d the return mapping algorithm which consists of pred ictor and corrector steps. The Cauchy stress at 𝑡௡ is expressed by the sum of Eq. (6) and Eq. (7) as follows, 𝝈ሺ𝑡௡ሻൌ𝑪 ஶ:𝜺𝒗𝒆ሺ𝑡௡ሻ൅෍𝒔 𝒊ሺ𝑡௡ሻூ ௜ୀଵ൅෍𝜎 ு௝ሺ𝑡௡ሻ𝟏௃ ௝ୀଵ (27) 11 where 𝑪 ஶൌ2 𝐺 ஶ𝑰ௗ௘௩൅3 𝐾 ஶ𝑰௩௢௟. First, by considering the elas tic deformation only, the predict or stress, 𝝈𝒕𝒓𝒊𝒂𝒍, is obtained from Eq. (27) and Eqs. (16, 17) as follows, 𝝈𝒕𝒓𝒊𝒂𝒍ሺ𝑡௡ାଵሻൌ𝑪 ஶ:𝜺𝒗𝒆ሺ𝑡௡ሻ൅𝑪෩:∆ 𝜺൅෍e x p ൬െ∆𝑡 𝑔௜൰𝒔𝒊ሺ𝑡௡ሻூ ௜ୀଵ ൅෍e x pቆ െ∆𝑡 𝑘௝ቇ𝜎 ு௝ሺ𝑡௡ሻ௃ ௝ୀଵ𝟏 (28) where 𝑪 ෩ൌ2 𝐺 ෨𝑰ௗ௘௩൅3 𝐾 ෩𝑰௩௢௟. Then, one determines if the plastic deformation is involved by considering the trial yield function, 𝑓௧௥௜௔௟, as follows, 𝑓௧௥௜௔௟ൌ𝜎 ௘௤௧௥௜௔௟െቀ 𝜎 ௬ሺ𝜺ሶሻ൅𝑅ሺ𝑝ሻቁw h e r e 𝜎 ௘௤௧௥௜௔௟ൌඨ3 2𝝈𝒕𝒓𝒊𝒂𝒍:𝝈𝒕𝒓𝒊𝒂𝒍 (29) If 𝑓௧௥௜௔௟൏0, deformation remains elastic region ( ∆𝜺𝒗𝒑ൌ0) and the corrector step is not necessary ( 𝝈ሺ𝑡௡ାଵሻൌ𝝈𝒕𝒓𝒊𝒂𝒍ሺ𝑡௡ሻ). If 𝑓௧௥௜௔௟൐0, which implies plastic deformation is involved, Cauchy stress at 𝑡௡ାଵ is adjusted as the correction o f predictor stre ss as follows, 𝝈ሺ𝑡௡ାଵሻൌ𝝈𝒕𝒓𝒊𝒂𝒍ሺ𝑡௡ሻെ𝑪෩:∆ 𝜺𝒗𝒑. (30) Combining Eq. (27) and Eq.(30), stress increment ( ∆𝝈 ≡ 𝝈 ሺ𝑡௡ାଵሻെ𝝈ሺ𝑡௡ሻ) over time step ( 𝑡௡, 𝑡௡ାଵ) is expressed as follows, ∆𝝈 ൌ 𝑪 ෩:ሺ∆𝜺 െ ∆𝜺𝒗𝒑ሻെ෍൤ 1െe x p ൬െ∆𝑡 𝑔௜൰൨𝒔𝒊ሺ𝑡௡ሻூ ௜ୀଵ െ෍ ቈ1െe x pቆ െ∆𝑡 𝑘௝ቇ቉𝜎ு௝ሺ𝑡௡ሻ𝟏௃ ௝ୀଵ (31) For the evaluation of the stress increment ( ∆𝝈) in Eq. (31), we approximate the viscoplastic strain increment ( ∆𝜺𝒗𝒑) from the flow rule (Eq. (23)) as 12 ∆𝜺𝒗𝒑ൎ∆ 𝑝 𝑵 . (32) For the unknown ∆𝑝 and 𝜎௘௤, Eq. (33) and Eq. (34) should be sat isfied at the same time. ∆𝑝 െ 𝑔 ௩∆𝑡 ൌ 0 (33) 𝜎௘௤൅3 𝐺 ෨∆𝑝 െ 𝜎 ௘௤௧௥௜௔௟ൌ0 (34) The accumulated plastic strain increment, ∆𝑝, is now can be obtained by iteratively solving Eq. (33) and Eq. (34) with New ton’s method (Doghri, 2013). 3. Adaptive incrementally affine linearized homogenization In the present section, by improvi ng the original incrementally affine scheme, we aim to propose a physically justifiabl e homogenization scheme based on the algorithmic tangent operator ( 𝑪𝟎𝒂𝒍𝒈, 𝑪𝟏𝒂𝒍𝒈) without isotropization or regularization. In the incrementall y affine method, stress for each phase ( 𝝈𝟎, 𝝈𝟏) are updated by accumulating their increments ( ∆𝝈 𝟎, ∆𝝈 𝟏) at each time step. Stress increments ( ∆𝝈 𝟎, ∆𝝈 𝟏) are determined by tangent operators at 𝑡௡ାଵ (𝑪𝟎𝒂𝒍𝒈ሺ𝑡௡ାଵሻ, 𝑪𝟏𝒂𝒍𝒈ሺ𝑡௡ାଵሻ) and rather independent to the history of tangent operators (𝑪𝟎𝒂𝒍𝒈ሺ𝑡ሻ, 𝑪𝟏𝒕𝒂𝒍𝒈ሺ𝑡ሻ 𝑡∈ሺ0, 𝑡 ௡ሻ) . However, the inconsistency between accumulated stresses (𝝈𝟎ሺ𝑡௡ሻ, 𝝈𝟏ሺ𝑡௡ሻ) due to the change in tangent operators and the shape change o f reinforcing particles is not taken into acc ount in previous studies. We int roduce the adaptive strain to correct the inconsistency. 3.1. Adaptive incrementally affine method with rigorous derivat ion The previous study derived the relationship between stress incr ement ( ∆𝝈) and strain increment ( ∆𝜺) based on the incrementally affine linearization and defined t he affine strain increment ( ∆𝜺𝒂𝒇) as follows (Miled et al., 2011), 13 ∆𝝈 ൌ 𝑪𝒂𝒍𝒈ሺ𝑡௡ାଵሻ:൫ ∆ 𝜺െ∆ 𝜺𝒂𝒇൯. (35) Algorithmic tangent operator ( 𝑪𝒂𝒍𝒈) is expressed as follows, 𝑪𝒂𝒍𝒈ൌ𝑪෩െ൫2𝐺෨൯𝟐 ℎ௩𝑵⊗𝑵െ൫ 2 𝐺 ෨൯𝟐 𝜎௘௤∆𝑝 𝜎௘௤൅3 𝐺 ෨∆𝑝𝜕𝑵 𝜕𝝈െ2𝐺෨ ℎ௩𝑔,ఙ𝑵⊗𝑔 ,𝜺 (36) where ℎ ௩≡1 ሺ∆𝑡ሻ𝑔,ఙ൅3 𝐺 ෨െ𝑔,௣ 𝑔,ఙ (37) 𝑔,ఙ≡𝜕𝑔௩ 𝜕𝜎௘௤;𝑔 ,௣≡𝜕𝑔௩ 𝜕𝑝;𝑔 ,ఌ≡1 ∆𝑡𝜕𝑔௩ 𝜕𝜺ሶ; (38) When the viscoplastic function ( 𝑔௩) is given as Eq. (26), Eq. (38) is arranged as follows. 𝑔,ఙൌ𝑚𝑔௩ 𝑓; 𝑔 ,௣ൌെ 𝑚 𝑔 ௩𝑑𝑅 𝑑𝑝ቆ1 𝑓൅1 𝜎௬൅𝑅ሺ𝑝ሻቇ; ℎ௩ൌ𝑓 ሺ∆𝑡ሻ𝑚𝑔 ௩൅3 𝐺 ෨൅𝑑𝑅 𝑑𝑝𝜎௘௤ 𝜎௬൅ 𝑅ሺ𝑝ሻ; (39) Affine strain increment ( ∆𝜺𝒂𝒇) is expressed as follows, ∆𝜺𝒂𝒇≡∆ 𝜺𝐞𝐯𝐩𝒂𝒇൅𝑪෩ି𝟏:ቐ෍൤ 1െe x p ൬െ∆𝑡 𝑔௜൰൨𝒔𝒊ሺ𝑡௡ሻூ ௜ୀଵ ൅෍ ቈ1െe x pቆ െ∆𝑡 𝑘௝ቇ቉𝜎ு௝ሺ𝑡௡ሻ𝟏௃ ௝ୀଵቑ (40) Here, the first term of Eq. (40) is related to inelastic strain increment and the second term of Eq. (40) is related to viscoel astic relaxation. In classical 𝐽ଶ viscoplastic model, ∆𝜺𝐞𝐯𝐩𝒂𝒇 i s expressed as follows, ∆𝜺𝐞𝐯𝐩𝒂𝒇ൌ𝑔 ௩ሺ𝑡௡ሻ∆𝑡ቈ𝑵ሺ𝑡௡ሻ൅𝑵 ሺ𝑡௡ାଵሻ𝑔,௣ሺ𝑡௡ାଵሻ∆𝑡 1െ𝑔 ,௣ሺ𝑡௡ାଵሻ∆𝑡቉. (41) In incrementally affine formulation, linearized homogenization is developed usi ng the analogy between thermo-elastic and viscoelastic-viscoplastic constituti ve equations. Having obtained the linearized relationship, the Mori-Tanaka method, a represen tative mean-field homogenization scheme, is applied. When affine strain increment s (∆𝜺𝟎𝒂𝒇, ∆𝜺𝟏𝒂𝒇) are zero, local 14 concentration tensor, 𝑩ଵ𝜺, which links the averaged strain increment of the matrix to th e averaged strain increment of the reinforcing particles is expre ssed as follows, 〈∆𝛆〉୵భൌ𝑩ଵ𝜺:〈∆𝛆〉୵బwhen ∆𝜺𝟎𝒂𝒇ൌ∆ 𝜺𝟏𝒂𝒇ൌ𝟎 (42) 𝑩ଵ𝜺ൌ൤ 𝑰൅𝑺 : ቀ 𝑪𝟎𝒂𝒍𝒈ቁିଵ ∶ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁ൨ିଵ (43) where 〈•〉୵బ, 〈•〉୵భ refer to the volume average over the matrix and particles, res pectively. Subscript zero and one refer to the matrix and particles, respe ctively (ex. 𝑪𝟎𝒂𝒍𝒈 refers to the algorithmic tangent operator of the matrix and 𝑪𝟏𝒂𝒍𝒈 refers to the algorithmic tangent operator of the particle). 𝑺 is the Eshelby’s tensor which i s a function of material propert ies of the matrix and shape of the particle. The relationship between Eshe lby’s tensor ( 𝑺) and Hill’s polarization tensor ( 𝑷) is defined as follows (Masson, 2008), 𝑃௜௝௞௟≡ቀ 𝑺 : 𝑪𝟎𝒂𝒍𝒈ି𝟏 ቁ ௜௝௞௟ൌ1 4𝜋නන 𝑀 ௜௝௞௟sin 𝜙 d𝜃d𝜙ଶగ ଴గ ଴ (44) where 𝑀 ௜௝௞௟ൌ1 4൫𝑁௝௞ିଵ𝜂௜𝜂௟൅𝑁௜௞ିଵ𝜂௝𝜂௟൅𝑁௝௟ିଵ𝜂௜𝜂௞൅𝑁௜௟ିଵ𝜂௝𝜂௞൯ 𝜂ଵൌsin 𝜙 𝑐𝑜𝑠 𝜃 𝑎ଵ;𝜂 ଶൌsin 𝜙 sin 𝜃 𝑎ଶ;𝜂 ଷൌcos 𝜙 𝑎ଷ; 𝑁௜௞ൌቀ 𝑪𝟎𝒂𝒍𝒈ቁ ௜௝௞௟𝜂௝𝜂௟;𝑥ଵଶ 𝑎ଵଶ൅𝑥ଶଶ 𝑎ଶଶ൅𝑥ଷଶ 𝑎ଷଶൌ1 where 𝑎ଵ, 𝑎ଶ, and 𝑎ଷ are the length of the half-axis of the ellipsoidal particle. G lobal concentration, 𝑨ଵ𝜺, which links macro st rain increment ( ∆𝜺̅) to strain average increment of the particle with zero affine strai n increments is expressed as fol lows, 𝑨ଵ𝜺ൌ𝑩ଵ𝜺:ሺ𝑣଴𝑰൅𝑣 ଵ𝑩ଵ𝜺ሻିଵ (45) where 𝑣଴ and 𝑣ଵ are volume fraction of the matri x and particles. Strain increm ents of the matrix and particle satisf y the following equation, ∆𝜺തൌ𝑣 ଴〈∆𝛆〉୵బ൅𝑣 ଵ〈∆𝛆〉୵భ. (46) The idea of this study is to corr ect the inconsistency between accumulated stresses ( 𝝈𝟎ሺ𝑡௡ሻ, 15 𝝈𝟏ሺ𝑡௡ሻ) due to the change in tangent operators and the shape change o f reinforcing particles. Each phases has stress and strain state at 𝑡௡ as as shown in Fig. 1 (a). In the sense of linearized comparison composite (LCC), stress and strain state can be rela ted as follows, 𝝈𝒊ሺ𝑡௡ሻൌ𝑪𝒊𝒂𝒍𝒈:𝜺𝒊ሺ𝑡௡ሻെ𝝈𝒊𝒂𝒇ሺ𝑡௡ሻൌ𝑪𝒊𝒂𝒍𝒈:ቀ 𝜺 𝒊ሺ𝑡௡ሻെ𝜺𝒊𝒂𝒇ሺ𝑡௡ሻቁ (47) where 𝝈𝒊𝒂𝒇ሺ𝑡௡ሻൌ𝑪𝒊𝒂𝒍𝒈:𝜺𝒊𝒂𝒇ሺ𝑡௡ሻ. (48) Subscript 𝑖 can be zero or one, which refers to state variable for the mat rix and particle, respectively. We note that the sum of the affine strain increme nt (∆𝜺𝒊𝒂𝒇) is not equal to the affine strain at 𝑡௡ because the modulus changes with strain. 𝜺𝒊𝒂𝒇ሺ𝑡௡ሻ്෍∆ 𝜺𝒊𝒂𝒇ሺ𝑡௞ሻ𝒏 𝒌ୀ𝟏 (49) Strain increment ( ∆𝜺𝒊) of each phase can be decomposed by strain increment in increm entally affine method ( ∆𝜺𝒊𝒊𝒂) and strain increment in adaptive scheme ( ∆𝜺𝒊𝒂𝒅) as follows, ∆𝜺𝒊ൌ𝜺 𝒊ሺ𝑡௡ାଵሻെ𝜺 𝒊ሺ𝑡௡ሻൌ∆ 𝜺𝒊𝒊𝒂൅∆ 𝜺𝒊𝒂𝒅. (50) ∆𝜺𝒊𝒊𝒂 is the term related to macro strain increment ( ∆𝜺ത), and ∆𝜺𝒊𝒂𝒅 is the term related to correction of inconsistency between accumulated stresses. Simil arly, the stress increment ( ∆𝝈 𝒊) of each phase can be decomposed by stress increment in incremen tally affine method ( ∆𝝈𝒊𝒊𝒂) and strain increment in adaptive scheme ( ∆𝝈𝒊𝒂𝒅) as follows, ∆𝝈 𝒊ൌ𝝈 𝒊ሺ𝑡௡ାଵሻെ𝝈 𝒊ሺ𝑡௡ሻൌ∆ 𝝈𝒊𝒊𝒂൅∆ 𝝈𝒊𝒂𝒅. (51) Strain increment of particle in incrementally affine scheme, ∆𝜺𝟏𝒊𝒂, is expressed as follows, and detailed explanation can be f ound in previous study (Miled et a l., 2011). ∆𝜺𝟏𝒊𝒂ൌ𝑨ଵ𝜺:∆ 𝜺ത൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ ∆ 𝝈𝟏𝒊𝒂െ𝑪𝟏𝒂𝒍𝒈:∆ 𝜺𝟏𝒊𝒂െ∆ 𝝈𝟎𝒊𝒂 ൅𝑪𝟎𝒂𝒍𝒈:∆ 𝜺𝟎𝒊𝒂ቁ (52) In adaptive scheme, we aim to obtain strain which satisfy the s train concentration with updated 16 algorithmic tangent operators and shape of particles. Consideri ng stress-strain state of Eq. (47), the strain concentration is calculated by following four steps, using the mathematical analogy between thermo-elastic and vi scoelastic-viscoplastic. Step 1. Assuming affine stresses of each phase ( 𝝈𝒊𝒂𝒇ሺ𝑡௡ሻ) are zero and macro strain 𝜺തሺ𝑡௡ሻ is subjected to composite, strain of the particle is obtained as follows, 𝜺𝟏𝒔𝟏ሺ𝑡௡ሻൌ𝑨ଵ𝜺:𝜺തሺ𝑡௡ሻ. (53) Step 2. In the presence of affine stresses, strain increment, ∆𝜺𝒔𝟐, is imposed to each phase so that each phase has uniform str ain and uniform stress state as follows, ∆𝝈 ൌ 𝑪𝟎𝒂𝒍𝒈:∆ 𝜺𝒔𝟐െ𝝈𝟎𝒂𝒇ሺ𝑡௡ሻൌ𝑪𝟏𝒂𝒍𝒈:∆ 𝜺𝒔𝟐െ𝝈𝟏𝒂𝒇ሺ𝑡௡ሻ (54) ∆𝜺𝒔𝟐ൌቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁି𝟏 :ቀ 𝝈𝟏𝒂𝒇ሺ𝑡௡ሻെ𝝈𝟎𝒂𝒇ሺ𝑡௡ሻቁ. (55) Step 3. Composite is subjected to macro strain increment, ∆𝜺𝒔𝟑, to eliminate the macro strain increment of step 2, which implies ∆𝜺𝒔𝟑ൌെ ∆ 𝜺𝒔𝟐. ∆𝜺𝟏𝒔𝟑ൌ𝑨ଵ𝜺:∆ 𝜺𝒔𝟑ൌെ 𝑨ଵ𝜺:∆ 𝜺𝒔𝟐 (56) Step 4. Superposing the strains of step 1-3, strai n of particle satisfies followi ng relationship. 𝜺𝟏ሺ𝑡௡ሻൌ𝑨ଵ𝜺:𝜺തሺ𝑡௡ሻ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ሺ𝑪𝟏െ𝑪 𝟎ሻି𝟏:ቀ 𝝈 𝟏ሺ𝑡௡ሻെ𝑪𝟏𝒂𝒍𝒈:𝜺𝟏ሺ𝑡௡ሻെ𝝈 𝟎ሺ𝑡௡ሻ ൅𝑪𝟎𝒂𝒍𝒈:𝜺𝟎ሺ𝑡௡ሻቁ (57) Here, algorithmic tangent operators ( 𝑪𝟎𝒂𝒍𝒈, 𝑪𝟏𝒂𝒍𝒈) in Eqs. (53-57) are defined in time step (𝑡௡ିଵ,𝑡௡). On the other hand, the algorithmic tangent operators and sha pe of particle are updated in time step ( 𝑡௡,𝑡௡ାଵ). The idea of this study is to correct the inconsistency betwe en accumulated stresses ( 𝝈𝟎ሺ𝑡௡ሻ, 𝝈𝟏ሺ𝑡௡ሻ) by considering the strain adj ustments in the matrix and particle ( ∆𝜺𝟎𝒂𝒅,∆ 𝜺𝟏𝒂𝒅). The unknown strain adjustments ( ∆𝜺𝟎𝒂𝒅,∆ 𝜺𝟏𝒂𝒅) satisfiy following relationship with updated algorith mic tangent operator and shap e of particle. 17 𝜺𝟏ሺ𝑡௡ሻ൅∆ 𝜺𝟏𝒂𝒅ൌ𝑨ଵ𝜺:𝜺തሺ𝑡௡ሻ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ሻ൅∆ 𝝈𝟏𝒂𝒅 െ𝑪𝟏𝒂𝒍𝒈:൫ 𝜺 𝟏ሺ𝑡௡ሻ൅∆ 𝜺𝟏𝒂𝒅൯െ𝝈 𝟎ሺ𝑡௡ሻെ∆ 𝝈 𝟎൅𝑪𝟎𝒂𝒍𝒈:൫ 𝜺 𝟎ሺ𝑡௡ሻ൅∆ 𝜺𝟎𝒂𝒅൯ቁ (58) where ∆𝝈𝟎𝒂𝒅 and ∆𝝈𝟏𝒂𝒅 are stress changes in adaptive scheme. Here, the strain changes for each phase satisfy following condition so that the macro strain increment in adaptive scheme is zero. ∆𝜺ത𝒂𝒅ൌ𝑣 ଴∆𝜺𝟎𝒂𝒅൅𝑣 ଵ∆𝜺𝟏𝒂𝒅ൌ0 (59) Combining Eq. (52) and Eq. (58), strain of the particle at 𝑡௡ାଵ (𝜺𝟏ሺ𝑡௡ାଵሻ) is obtained as follows, 𝜺𝟏ሺ𝑡௡ାଵሻൌ𝜺 𝟏ሺ𝑡௡ሻ൅∆ 𝜺 𝟏 ൌ𝑨ଵ𝜺:𝜺തሺ𝑡௡ାଵሻ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ାଵሻെ𝑪𝟏𝒂𝒍𝒈:𝜺𝟏ሺ𝑡௡ାଵሻെ𝝈 𝟎ሺ𝑡௡ାଵሻ ൅𝑪𝟎𝒂𝒍𝒈:𝜺𝟎ሺ𝑡௡ାଵሻቁ. (60) Eq. (60) is an implicit equation of ∆𝜺𝟏 because strain incrment of the matrix is determined by Eq. (46) and the algorithmic stif fness tensors and stress incre ments for each phase are determined by its strain increments. So, ∆𝜺𝟏 is obtained by backward iterative method and then other variables are determined. Macro stress increment is computed using the strain increments of each phase as follows, ∆𝝈ഥൌ𝑣 ଴𝑪𝟎𝒂𝒍𝒈:ቀ〈∆𝜺〉௪బെ∆ 𝜺𝟎𝒂𝒇ቁ൅𝑣 ଵ𝑪𝟏𝒂𝒍𝒈:ቀ〈∆𝜺〉௪భെ∆ 𝜺𝟏𝒂𝒇ቁ. (61) The effective algorithmic ta ngent operator of composite, 𝑪ഥ, is expressed as follows, 𝑪ഥൌቂ 𝑣 ଴𝑪𝟎𝒂𝒍𝒈൅𝑣 ଵ𝑪𝟏𝒂𝒍𝒈:𝑩ଵ𝜺ቃ:ሾ𝑣଴𝑰൅𝑣 ଵ𝑩ଵ𝜺ሿିଵ. (62) To summarize, the adaptive incrementally affine homogenization proceeds following orders: 1. Initialize the strain increme nt of the particles as macro st rain increment 〈∆𝜺〉୵భൌ∆ 𝜺 ̅ 2. Compute the algorithmic tangent operator ( 𝑪𝟏𝒂𝒍𝒈) from Eq. (36) using 〈∆𝜺〉୵భ and 18 〈𝜺ሺ𝑡 ௡ሻ〉௪భ as inputs. 3. The strain increment of the m atrix is calculated from Eq. (4 6) 4. Compute the algorithmic tangent operator ( 𝑪𝟎𝒂𝒍𝒈) from Eq. (36) using 〈∆𝜺〉௪బ and 〈𝜺ሺ𝑡 ௡ሻ〉௪బ as inputs. 5. Eshelby’s tensor ( 𝑺) and local concentration tensor ( 𝑩ଵ𝜺) are computed from Eq. (44) and (43) using 𝑪𝟎𝒂𝒍𝒈. 6. Residual tensor 𝑹 is computed as follows, 𝑹ൌ 𝑨ଵ𝜺:𝜺തሺ𝑡௡ାଵሻെ𝜺 𝟏ሺ𝑡௡ାଵሻ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ାଵሻെ𝑪𝟏𝒂𝒍𝒈:𝜺𝟏ሺ𝑡௡ାଵሻെ𝝈 𝟎ሺ𝑡௡ାଵሻ ൅𝑪𝟎𝒂𝒍𝒈:𝜺𝟎ሺ𝑡௡ାଵሻቁ. (63) 7-1. If |𝑹|൐t o l e r a n c e , go to step 2 with a new trial from Eq. (64). 〈𝜺𝟏ሺ𝑡௡ାଵሻ〉௪భ←〈𝜺𝟏ሺ𝑡௡ାଵሻ〉௪భ൅𝑹 (64) 7-2. If |𝑹|൏t o l e r a n c e , then go to step 8. 8. Macro stress increment and effective tangent operator are ob tained from Eq. (61) and Eq. (62) 3.2. Adaptive incrementally affine method for better prediction The adaptive incrementally affine method introduced in section 3.1 is mathematically rigorous, but it seems to involve some errors in the plastic deformation regime with finite particle volume fraction. In this section, we woul d like to sugg est a method fo r better predict ion by making an assumption. As shown in Fig. 1 (a), linearized state includes a ffine strain and affine stress. In this section, when we take the li nearization for each phase, we made an assumption excluding affine strain and affine stress a s s h o w n i n F i g . 1 ( b ) . A t t h i s t i m e , t h e a l g o r i t h m i c t a n g e n t 19 operator can be seen as a secant tangent operator considering t he newly defined strain 𝑪𝒊𝒂𝒍𝒈ି𝟏 :𝝈𝒊ሺ𝑡௡ሻ. With the assumption, adaptive strain increments ( ∆𝜺𝟎𝒂𝒅, ∆𝜺𝟏𝒂𝒅) are derived in the following procedures. Step 1. When a fictitious strain increment ∆𝜺𝒔𝟏 is added to each phase at 𝑡௡ with updated 𝑪𝒂𝒍𝒈ሺ𝑡௡ାଵሻ, fictitious stresses ( 𝝈𝟎∗, 𝝈𝟏∗) are evolved as follows, 𝝈𝟎∗ൌ𝑪𝟎𝒂𝒍𝒈ሺ𝑡௡ାଵሻ:∆ 𝜺𝒔𝟏൅𝝈 𝟎ሺ𝑡௡ሻ (67) 𝝈𝟏∗ൌ𝑪𝟏𝒂𝒍𝒈ሺ𝑡௡ାଵሻ:∆ 𝜺𝒔𝟏൅𝝈 𝟏ሺ𝑡௡ሻ (68) At this time, ∆𝜺𝒔𝟏 is set to satisfy the followi ng local concentration condition, 𝑪𝟏𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟏∗ൌ𝑩ଵ𝜺:𝑪𝟎𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟎∗ (69) Substituting the Eq. (67) and Eq. (68) into Eq. (69), the ficti tious strain ( ∆𝜺𝒔𝟏) is expressed as follows, ∆𝜺𝒔𝟏ൌሺ𝑰െ𝑩ଵ𝜺ሻିଵ:ቂ 𝑩ଵ𝜺:𝑪𝟎𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟎ሺ𝑡௡ሻെ𝑪𝟏𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟏ሺ𝑡௡ሻቃ (70) Step 2. ∆𝜺𝒔𝟐ሺൌ െ∆𝜺𝒔𝟏ሻ is applied to the composite seque ntially, so that macro strain increment by fictitious strain is zero. When ∆𝜺𝒔𝟐 is applied to the composite, 𝑨ଵ𝜺:∆ 𝜺𝒔𝟐 o f strain change occurs in particl es. The procedure of step 1 to s tep 2 is illustrated in Fig. 2. ∆𝜺𝟏𝒔𝟐ൌ𝑨ଵ𝜺:∆ 𝜺𝒔𝟐ൌെ 𝑨ଵ𝜺:∆ 𝜺𝒔𝟏 (71) Step 3. Considering Eq. (55) and Eq. (56), strain increment of the particle by affine strain increment in adaptive sche me is obtained as follows, ∆𝜺𝟏𝒔𝟑ൌሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ሻ൅∆ 𝝈𝟏𝒂𝒅െ𝑪𝟏𝒂𝒍𝒈:൫ 𝜺 𝟏ሺ𝑡௡ሻ൅∆ 𝜺𝟏𝒂𝒅൯ െ𝝈 𝟎ሺ𝑡௡ሻെ∆ 𝝈 𝟎൅𝑪𝟎𝒂𝒍𝒈:൫ 𝜺 𝟎ሺ𝑡௡ሻ൅∆ 𝜺𝟎𝒂𝒅൯ቁ (72) Step 4. Strain incremental of the particle in adaptive scheme i s obtained by superposing step 1-3 as follows, 20 ∆𝜺𝟏𝒂𝒅ൌሺ𝑰െ𝑨ଵ𝜺ሻ:ሺ𝑰െ𝑩ଵ𝜺ሻିଵ:ቂ 𝑩ଵ𝜺:𝑪𝟎𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟎ሺ𝑡௡ሻെ𝑪𝟏𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟏ሺ𝑡௡ሻቃ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ሻ൅∆ 𝝈𝟏𝒂𝒅 െ𝑪𝟏𝒂𝒍𝒈:൫ 𝜺 𝟏ሺ𝑡௡ሻ൅∆ 𝜺𝟏𝒂𝒅൯െ𝝈 𝟎ሺ𝑡௡ሻെ∆ 𝝈 𝟎൅𝑪𝟎𝒂𝒍𝒈:൫ 𝜺 𝟎ሺ𝑡௡ሻ൅∆ 𝜺𝟎𝒂𝒅൯ቁ (73) Then, combining Eq. (52) and Eq. (73), strain increment of the particle during time step (𝑡௡,𝑡௡ାଵ) is expressed as follows, 〈∆𝜺〉௪భൌ𝑨ଵ𝜺:∆ 𝜺ത ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ାଵሻെ𝑪𝟏𝒂𝒍𝒈:𝜺𝟏ሺ𝑡௡ାଵሻെ𝝈 𝟎ሺ𝑡௡ାଵሻ ൅𝑪𝟎𝒂𝒍𝒈:𝜺𝟎ሺ𝑡௡ାଵሻቁ ൅ሺ𝑰െ𝑨ଵ𝜺ሻ:ሺ𝑰െ𝑩ଵ𝜺ሻିଵ:ቂ 𝑩ଵ𝜺:𝑪𝟎𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟎ሺ𝑡௡ሻ െ𝑪𝟏𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟏ሺ𝑡௡ሻቃ (74) Here, only difference from Section 3.1 is that Eq. (75) and Eq. (76) are used instead of Eq. (63) and Eq. (64) for the proce dure explained in Section 3.1. 𝑹ൌ 𝑨ଵ𝜺:∆ 𝜺തെ〈∆𝜺〉௪భ ൅ሺ𝑨ଵ𝜺െ𝑰ሻ:ቀ 𝑪𝟏𝒂𝒍𝒈െ𝑪𝟎𝒂𝒍𝒈ቁିଵ :ቀ 𝝈 𝟏ሺ𝑡௡ାଵሻെ𝑪𝟏𝒂𝒍𝒈:𝜺𝟏ሺ𝑡௡ାଵሻെ𝝈 𝟎ሺ𝑡௡ାଵሻ ൅𝑪𝟎𝒂𝒍𝒈:𝜺𝟎ሺ𝑡௡ାଵሻቁ ൅ሺ𝑰െ𝑨ଵ𝜺ሻ:ሺ𝑰െ𝑩ଵ𝜺ሻିଵ:ቂ 𝑩ଵ𝜺:𝑪𝟎𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟎ሺ𝑡௡ሻ െ𝑪𝟏𝒂𝒍𝒈ିଵሺ𝑡௡ାଵሻ:𝝈𝟏ሺ𝑡௡ሻቃ (75) 〈∆𝜺〉௪భ←〈∆𝜺〉௪భ൅𝑹 (76) Adaptive incrementally affine method explained in this section is used for verification of Section 4 instead of mathematica lly rigorous formulas in Sectio n 3.1. 3.3. Yield reduction method Despite the elaborate formulation, we still find that the plast ic deformation in FEM simulations is initiated earlier than in homogenization based on the adapti ve incrementally affine method. This is expected because, while non-uniform stress occurs withi n the matrix around the particle-matrix interface in FEM, the strain/stress is determin ed from the average strain/stress 21 over the entire matrix region in mean-field homogenization. Hen ce, in the homogenization scheme, plastic deformation is not considered until the average d mean-field stress satisfies the yield condition, while the region with local stress concentrati on undergoes the plastic deformation in FEM simulations. The non-uniformity of stress-fi eld increases with the particle volume fraction due to enhanced particle-particle interaction. Hence, the greater discrepancy exists between the homogeniza tion and FEM results for composite s involving more reinforcing particles. To resolve the probl em, we propose the yield reducti on method (Eq. 77) which enforce the earlier initiation of the plastic deformation in th e homogenization scheme as follows, 𝜎୷←ሺ1െ𝑣 ଵሻ𝜎୷ for the matrix in homogenization (77) A greater reduction of yield str ess is applied with a larger vo lume fraction of particles to account for enhanced particle-par ticle interaction. The yield r eduction method is applied only to the matrix not to the particles because the stress/strain fi eld within ellipsoidal/spherical inclusions are almost uniform under the elastic regime and even after the initiation of the local plastic deformation 4. Result and discussion We test the validity of our homogenization scheme against the F EM simulations of particle- reinforced composites under comp lex loading conditions by varyi ng the volume fraction and shape of particles. The commercial software ABAQUS (Abaqus, 201 7) with UMAT code is used for FEM. FCC structured unit cells which have the ideal di spersion of particles in three- dimensional space are used as shown in Fig. 3. Periodic boundar y conditions are applied to every face of unit cells. For the unit cell with spherical part icles (Fig. 3 (a)), 220,779 and 22 243,717 of quadratic tetrahedral elements are used for 5% and 1 5% of volume fraction, respectively. For the unit cell with ellipsoidal particles (Fig . 3 (b)), 343,311 of quadratic tetrahedral elements are used for 15% of volume fraction. Our h omogenization scheme is implemented as the UMAT code in ABAQUS. Strain, strain incremen t, and stress are passed to UMAT at every step, and then UMAT returns the Cauchy stress and tangent operator of the next step according to the adaptive incremental affine method. While FEM is performed with numerous elements, a single brick element is used to obtain the prediction from the adaptive incrementally affine method. Various simulations are carried ou t for the cases considering viscoelastic, elastic-viscoplastic, and viscoelastic-viscoplast ic matrices under uniaxial, cyclic, and bi-axial loading. Different strain rate are used for each t ests (5 [1/s], 0.5 [1/s], 0.005 [1/s]) to test if the rate-dependency is well accounted for via our ho mogenization scheme. We note that verification is conducted u sing adaptive incrementally aff ine method introduced in Section 3.2. 4.1. Viscoelastic-viscoplastic matrix with elastic particles The material properties of the m atrix and particles are shown i n Table. 1. For the validation of viscoelastic behavior, a composite of the viscoelastic matri x reinforced with elastic spherical particles is tested. The yield stress of the matrix is set to i nfinity ( 𝜎௬→∞ ) to model viscoelastic behavior. Predictions from the adaptive incrementally affine me thod are compared to FEM results at 15% of the volume fraction in Fig. 4. While the incr ementally affine method showed stiff prediction in viscoelastic composite (Miled et al., 2013) , the adaptive incrementally affine method in Fig. 4 shows a very good match with FEM results (less than 1.3% of error). These results confirm that adaptive strain compensates for the incons istency of accumulated stress according to the shape change of the particles, which was not c onsidered in the previous 23 incrementally affine method. A composite of the elastic-visco plastic matrix reinforced with elastic spherical particles is tested. Relaxation moduli in Table. 1 are set to zero to model elastic-viscoplastic behavior ( 𝐺௜ൌ 𝐾௝ൌ0). Predictions are compared to F EM results at 5% and 15% of vol ume fractions in Fig. 5. Predictions show a good match at 5% of volume fraction, but a notable discrepancy between predictions and FEM is observed at 15% of volume fraction in th e plastic deformation region. On the other hand, the incrementally affine method with isotrop ization and regularization technique from the previous paper has some error even in small volume fraction (Fig. S.1). In the sense of mean-field homoge nization, predictions should have good accuracy at low volume fraction, and error increases with volume fraction due to the i nteraction between particles. Comparing Fig. 5 and Fig. S.1, the adaptive incrementally affin e method is assured to show better prediction in both small and finite volume fraction. Fig . 5 shows that the discrepancy between prediction and FEM results is caused because the plasti c deformation initiates earlier in FEM than in the homogenization. As shown in Fig. 6, in FEM, plastic deformation initiate in local areas due to stress concentration while most areas sti ll undergo elastic deformation. On the other hand, in homogenization, plastic deformation starts w hen mean-stress satisfies the yield condition. To compensate for the difference between FEM a nd homogenization, the yield reduction method (Eq. ( 77)) is applied to yi eld stress of the m atrix in homogeni zation. Results from the adaptive incrementally affine method with yield reduct ion method for different strain rates are shown in Fig. 7. Predictions and FEM results with a d ifferent viscoplastic exponent (𝑚) are compared in Fig. 8. Predictions show a very good match wi th FEM results (less than 2% of error). A maximum error is observed near the yield point. A composite of the viscoelastic-viscoplastic matrix reinforced with elastic spherical particles is tested. Fig. 9 (a, b) show predictions of the adaptive incre mentally affine method without 24 yield reduction method, and Fig. 9 (c, d) show predictions of t he adaptive incrementally affine method with yield reduction method. As we discussed before, pre dictions without yield reduction show better predictions at low volume fraction and no table errors at high volume fraction. On the other hand, pred ictions with yield reduction s how a good match with FEM results (less than 2.3% of error). Results for 20% and 30% of v olume fractions are shown in Fig. S.2 showing good prediction up to 20% of volume fraction. Results for adaptive incrementally affine with rigorous derivation explained in Sect ion 3.1 are shown in Fig. S.3. The adaptive incrementally affine with rigorous derivation seem s to involve error in plastic deformation regime. Relaxation tests are conducted fo r the viscoelastic-viscoplast ic matrix reinforced with elastic spherical particles as s hown in Fig. 10. The composite specimen is loaded up to 0.03 of strain with different strain rates and then relaxed. Uniaxial c yclic tests are performed for the viscoelastic-viscoplastic matr ix reinforced with elastic spheri cal particles as shown in Fig. 11. Even though the adaptive incrementally affine method without yi eld reduction shows fine results as shown in Fig. 11 (a), the adaptive incrementally aff ine method with yield reduction shows better predictions as shown in Fig. 11 (b). Predictions u nder two-cycle and three-cycle conditions are presented in Fig. S.4. Bi-axial tests are perfor med with different strain rates for the viscoelastic-viscoplastic matrix reinforced with elastic sp herical particles as shown in Fig. 12, Fig. 13, and Fig. 14. Strain rates in x-direction and y-direction are equal in Fig. 12. Strain rates in x-direction and y-direction have the same magnitude but opposite signs in Fig. 1 3. Strain rates in the y-direction are half of the strain rates of x-direction in Fig. 14. The adaptive incrementally affine method with yield reduction shows very goo d prediction under various loading conditions. Uniaxial loading test s are conducted for th e viscoelastic-visco plastic matrix reinforced with 25 elastic ellipsoidal particles as shown in Fig. 15. The aspect r atio of ellipsoidal particles is set to three and the unit cell in Fig. 3 (b) is used for FEM. Fig. 15 (a) shows the results when the x-directional loading is subjected to unit cell and Fig. 15 (b) shows the results when the y- directional loading is subjected to the unit cell. It is observ ed that predictions with ellipsoidal particles have some errors compared to the predictions with sph erical particles. Stresses of each phase are presented in Fig. 16 to investigate the origin of the errors. When the x-directional loading is subjected to the composite, prediction for the parti cle is a little soft. On the other hand, when the y-directional loading is subjected to the composite, predictions for the particle is a little stiff, and predicti on for the matrix is a little so ft. Notable errors are observed beyond the yield point. Violation of uniform field within the ellipsoi dal particles is suspected as a cause of the errors which infringe the assumption of the Mori-Tanaka method. However, the adaptive incrementally affine method still demonstrates much improved pr ediction compared to the incrementally affine method from the previous paper. 4.2. Elastic-viscoplastic matrix w ith elastic-viscoplastic part icles Uniaxial tests are conducted when both the matrix and particles have elastic-viscoplastic behavior. Material pr operties for the matrix and particles are shown in Table. 2. Fig. 17 shows the results for predictions with and without yield reduction. B oth methods show good prediction compared to FEM results (less than 2.3% of errors). Predictions with yield reduction are very close to the predicti ons without yield reduction (less than 0.6% of difference) because the ratio of yield stress to young’s modulus of the matrix is v ery small ( 𝜎௬/𝐸଴ൌ ሺ70 MPaሻ/ሺ70 GPaሻ ). Even so, predictions with the yield reduction method show cl oser prediction to FEM results. 26 5. Conclusion In this study, we have developed the adaptive incrementally aff ine method which is applicable for predicting the mechanical response of viscoelast ic-viscoplastic composite. We proposed an improved formulation of the incrementally affine me thod by considering the consistency of the concentrati on tensor and the stress state at each step of loading increment, which enables the correct prediction of the viscoelastic behavi or. However the rigorous formulas seems to involve some errors in plastic deformation re gime. So, we made an assumption dropping out the affine strain and affine stress in adaptive scheme, which leads to better prediction. Because the mean -field homogenization method does not fully account for the non-uniform distribution of stress field around the reinfor cements, the initiation of the plastic deformation in realistic composite occurs earlier than the homogenization method. Hence, in order to resolve the i ssue, we propose the yield redu ction method in which the yield stress of the matrix is assume d to be monotonically decreasing function of the reinforcement volume fraction. The prediction fr om the adaptive incremental a ffine homogenization scheme with yield reduction turns out to match well with the FEM simul ation of viscoelastic, elastic- viscoplastic, and viscoelastic-v iscoplastic matrix composites u nder various load ing conditions such as uniaxial, cycli c, and bi-axial loadings. 27