Despite intense research on the homogenization methods, it stil l is a challenging task to predict the nonlinear mechanical responses of visco-hyperelastic partic ulate-reinforced composites. In this work, we propose the adap tive affine met hod, a novel mean- field homogenization method designed to ensure the consistency of the accumulated strain st ate and the concentration tensor, and apply the method to predict the mechanical response of the composite in the large strain regime under uniaxial, cyclic, and bi-axial loadings. Our metho d is also extended to predict the mechanical response in the presen ce of interfacial imperfection s described by linear cohesive traction-separation laws. The ana lytic predictions are validate d against finite element analyses of representative volume elements. We believe that our adaptive affine method can be extended to model various nonlinear responses of load-bearing composites including the effects of (visco)plasticity and finite deformation. 3 1. Introduction Particle reinforced composites have been widely adopted to enha nce various physical properties including the structura l, thermal, thermoelectric, p iezoelectric properties, as well as to gain multi-functionality [1-4]. For the industrial developme nt of composites, it is a prerequisite to predict their properties in terms of the compos ition of their constituents and applied external mechanical, thermal, and electrical stimuli. I n that regard, the homogenization theory is well established for the prediction of effective phys ical properties of composites under a linear response regime. However, there is still room for impr ovement when it comes to the nonlinear response originated fro m geometric or material nonlin earity. One typical example of particle reinforced composite in which both geometric and mater ial nonlinearities play a critical role is elastomeric matrix composites. For example, solid prope llant materials, which consist of rubbery matrix and oxidizer pa rticles such as aluminum or am monium perchlorate, have been extensively studied due to their importance in designing h igh precision rocket and missile propulsion systems [5-10]. Besid es solid propellants, visco-hyp erelastic elastomer composites are widely used in stretchable electronics, soft robotics, wear able systems, automobile industries, etc. Accordingly, a v ariety of homogenization frame works that can appropriately account for such nonlinear res ponses have been suggested. Many previous researchers have adopted mean-field homogenizatio n (MFH) approaches to predict the effec tive properties of solid propell ants [11-16]. For the linear elastic response, one can apply the well-established MFH approaches suc h as the Mori-Tanaka (M-T) scheme which provides a closed-form solution for the effective elastic properties of the composites by utilizing Eshelby’s solution for the single inclu sion problem. However, for the prediction of inelastic behaviors, multiple methods were sugges ted to apply MFH and it remains controversial which frameworks provide numerically effi cient and mathematically rigorous predictions. A frequenc y-domain homogenization method was established for linear viscoelastic composites in which the correspondence principle h olds by utilizing the Laplace- Carson transform (LCT) [17]. However, the reconstruction of the solution in the time domain involves a costly numerical inversion of the LCT. A similar pro cedure is followed for the MFH of elasto-viscoplastic composites with an integral affine formu lation [18]. Besides being costly, the numerical inversion of the LCT by a collocation method migh t be inaccurate. Possible improvements have been studi ed by Rekik and Brenner [19]. 4 On the other hand, significant efforts have been devoted to dev eloping time-domain homogenization methods for the sake of mathematical simplicity and convenience. The linear comparison composite (LCC) concept was proposed to extend linea r MFH theories to nonlinear constituents by approaching the actual response within a given time step via virtual linear materials. The LCC can be defined based on variational formulat ions [20, 21] or the linearization of the local cons titutive models. Depending on th e linearization scheme, secant or tangent operators are adopted [22-24]. The incrementally affine method was developed using the analogy between linearized ela sto-viscoplasticity and therm oelasticity [25] and between linearized viscoelasticity-viscoplasticity and thermoelasticity [26]. However, the incrementally affine method tends to show stiff predictions in the viscoelast ic regime. The prediction becomes stiffer in visco-hyperelastic mat erial under finite deformation . Alternatively, the incremental secant method was developed whic h combines the advantages of secant modulus and tangent modulus methods [27-30 ]. The method showed acceptable predictability for mat erials with elasto-plasticity or elasto-viscoplasticity, but its extension to viscoelastic-visc oplastic composites has not been proposed yet. In this work, we propose a new homogenization scheme, named the adaptive affine method, that can resolve the aforementioned problems. We found that existing methods suffer from the inconsistency of stress state at 𝑡௡ when performing the homogenization of linearized material in a given time step [ 𝑡௡, 𝑡௡ାଵ], i.e. the stress concentration tensor does not correctly indicate the relationship betw een the applied stress and the st ress at the interior of reinforcements. To resolve the problem, we propose a method whi ch adaptively adjusts the strain of each constituent at eve ry step of the loading process to ensure the consistency of the accumulated strain state and the concentration tensor. We utili ze the 4th order tangent operator as it is, which enables us to analyze general loading cases inc luding cyclic loading and multi- axial loading under finite deformation. Our method is then exte nded to predict the mechanical response in the presence of the i nterfacial imperfections descr ibed by a linear cohesive traction- separation law. For the linear interfacial damage model, the in terfacial damage is represented as finite compliance of interf acial spring, and the homogenizat ion scheme is implemented using an effective inclusion method. Our homogenization model shows a good match with three- dimensional (3D) finite element model (FEM) simulations. 5 This paper consists of the following sections. The modeling of visco-hyperelastic material is reviewed based on Simo’s model in Section 2. The ad aptive affine method is introduced and compared with the incrementally affine method in Section 3. The adaptive affine method for interfacial damage cases is introduced in Sec tion 4. The conclusion is reached in Section 5. Boldface letters indicate second- or fourth-order tensors and n on-bold face letters indicate scalar values. (:) ref ers to double contraction and ( ⊗) refers to the dyadic product. 𝒂: 𝒃 ൌ 𝑎 ௜௝𝑏௝௜; ሺ𝑨:𝒃ሻ௜௝ൌ𝐴௜௝௞௟𝑏௟௞ሺ𝒂⊗𝒃 ሻ௜௝௞௟ൌ𝑎 ௜௝𝑏௞௟ 2. Visco-hyperelastic material modeling In this paper, we adopt Simo’s visco-hyperelastic model because of the following advantages [31-33]. First, the vol umetric or deviatoric part ca n be selectively chosen as a source of viscosity because the formulation separates the volum etric and deviatoric parts. Second, the model uses a form of strain energy density for the constitutive equation, which allows most of the existing hyperelastic model to be combined. Third, the model can be extended to reflect thermal eff ects by using the time-temperatu re superposition principle. Elastic stress is derived from strain energy density ( 𝜓∘) which is additively decomposed into volumetric part ( 𝑈∘) and deviatoric part ( 𝑊ഥ∘) as follows, 𝜓∘ൌ𝑈∘ሺ𝐽ሻ൅𝑊ഥ∘ሺ𝑪ഥሻ (1) Then, the second Piola-Kirchhoff stress ( 𝑺ሺ𝑡ሻ) is derived from strain energy density with a viscous response as follows, 𝑺ሺ𝑡ሻൌ𝐽 𝑃 𝑪ିଵ൅𝑯 (2) where 𝑃 ൌ න𝐾ሺ𝑡െ𝑠 ሻ 𝐾଴𝜕 𝜕𝑠൬𝜕𝑈∘ 𝜕𝐽൰𝑑𝑠௧ ିஶ, (3) 𝑯ൌන𝜇ሺ𝑡െ𝑠 ሻ 𝜇଴𝜕𝑺ഥ∘ 𝜕𝑠𝑑𝑠௧ ିஶ, (4) 𝑺ഥ∘ൌ𝐽ିଶ ଷ𝑺෨∘ൌ𝐽ିଶ ଷDEV ሺ2𝜕௖̅𝑊ഥ∘ሻ (5) 6 Here, 𝐾 and 𝜇 are bulk and shear moduli, respectively, 𝐽 is the determinant of deformation gradient ( 𝑭), and 𝑪 is the right Cauchy-Green strain tensor ( 𝑪) which is defined as 𝑪≡𝑭்𝑭. Deviatoric parts of the deformation gradient ( 𝑭ഥ) and right Cauchy-Green strain tensor ( 𝑪ഥ) are expressed as 𝑭ഥൌ𝐽ିଵ/ଷ𝑭 and 𝑪ഥ≡𝐽ିଶ/ଷ𝑪, respectively. 𝜕௖̅ refers to partial derivative with respect to 𝑪ഥ. DEV ሺ൉ሻ refers to deviatoric projection operator in the material descr iption which is defined as follows, DEV ሺ൉ሻൌሺ൉ሻെ1 3ሾሺ൉ሻ:𝑪ሿ𝑪ି𝟏 (6) Viscoelastic properties are re presented with the Prony series a s follows, 𝜇ሺ𝑡ሻൌ𝜇 ஶ൅෍𝜇 ௜expሺെ𝑡/𝜏 ௜ሻூ ௜ୀଵ (7) where 𝜇 ଴ൌ𝜇ሺ0ሻ, 𝑔ሺ𝑡ሻ≡ఓሺ௧ሻ ఓబ, 𝑔ஶ≡ఓಮ ఓబ, 𝑔௜≡ఓ೔ ఓబ 𝐾ሺ𝑡ሻൌ𝐾 ஶ൅෍𝐾 ௝exp൫െ𝑡/𝜂 ௝൯௃ ௝ୀଵ (8) where 𝐾 ଴ൌ𝐾 ሺ0ሻ In this paper, we take only the deviatoric part as a source of viscosity for simplicity, which implies that every 𝐾௝ is taken to be zero. The elastic tangent operator ( 𝑳∘) is expressed as follows, 𝑳∘ൌ𝑳∘ഥ൅4ቆ𝜕ଶ𝑈∘ 𝜕𝐽ଶ𝜕𝐽 𝜕𝑪⊗𝜕𝐽 𝜕𝑪൅𝜕𝑈∘ 𝜕𝐽𝜕ଶ𝐽 𝜕𝑪𝜕𝑪ቇ (9) 𝑳∘ഥൌ𝑳ௗ∘തതതെ2 3ቆ𝑺ഥ∘⊗𝑪ି𝟏൅𝑪ି𝟏⊗𝑺ഥ∘൅𝐽ିଶ/ଷሺ2𝜕௖̅𝑊ഥ∘:𝑪ሻ൬𝑰𝑪ష𝟏൅1 3𝑪ି𝟏⊗𝑪ି𝟏൰ቇ (10) 𝑳ௗ∘തതതൌ4 𝐽ିସ/ଷ൤𝜕𝑪ഥ𝑪ഥଶ𝑊ഥ∘൅1 9൫𝑪: 𝜕𝑪ഥ𝑪ഥଶ𝑊ഥ∘:𝑪 ൯ 𝑪ି𝟏⊗𝑪ି𝟏െ1 3𝑪ି𝟏⊗൫ 𝜕𝑪ഥ𝑪ഥଶ𝑊ഥ∘:𝑪 ൯ െ1 3൫𝜕𝑪ഥ𝑪ഥଶ𝑊ഥ∘:𝑪 ൯⊗𝑪ି𝟏൨ (11) 7 where 𝑰𝑪షభ≡𝜕ሺ𝑪ିଵሻூ௃ 𝜕𝑪 ௄௅ൌെ1 2ൣሺ𝑪ିଵሻூ௄ሺ𝑪ିଵሻ௃௅൅ሺ𝑪ିଵሻூ௅ሺ𝑪ିଵሻ௃௄൧ (12) The elastic tangent operator also can be considered as an insta ntaneous operator. The first term in Eq. (9) corresponds to the deviatoric part while the second term corresponds to the volumetric part. In the case of Neo-Hookean solids, 𝑳ௗ∘തതത goes to zero. Numerical calculation of the stress including the viscoelastic effect can be obtained as below. Substituting Eq. (7) into Eq. (4), second Piola-Kirchhoff stress at 𝑡௡ାଵ is expressed as follows, 𝑺௡ାଵ ൌ𝐽 ௡ାଵ𝜕𝑈∘ሺ𝐽௡ାଵሻ 𝜕𝐽𝑪௡ାଵିଵ൅𝑔 ஶ𝑺ഥ௡ାଵ∘ ൅෍𝑔 ௜ே ௜ୀଵ𝐽௡ାଵିଶ ଷDEV ௡ାଵሺ𝑯௡ାଵ௜ሻ (13) 𝑯௡ାଵ௜ is called an internal algorithmic variable which has an integr a l f o r m . T h e d e t a i l e d process to calculate 𝑯௜ is shown below. 𝑯௜ሺ𝑡ሻ≡න e x p ሾെሺ𝑡െ𝑠 ሻ/𝜏௜ሿ𝑑𝑺෨∘ 𝑑𝑠𝑑𝑠௧ ିஶ (14) 𝑯௡ାଵ௜ൌe x p ሺെΔ𝑡/𝜏 ௜ሻ𝑯௜ሺ𝑡௡ሻ൅න e x p௧೙శభ ௧೙ሾെሺ𝑡௡ାଵ െ𝑠ሻ/𝜏௜ሿ𝑑𝑺෨∘ 𝑑𝑠𝑑𝑠 (15) 𝑯௡ାଵ௜ൌe x p ሺെΔ𝑡/𝜏 ௜ሻ𝑯௡௜൅e x p ሺെΔ𝑡/2𝜏 ௜ሻ൫𝑺෨௡ାଵ∘െ𝑺෨௡∘൯ (16) 𝑯෩௜ is newly defined from 𝑯௜ to separate variables within the time step [ 𝑡௡, 𝑡௡ାଵ]. 𝑯෩௡௜≡e x p ሺെΔ𝑡/𝜏 ௜ሻ𝑯௡௜െe x p ሺെΔ𝑡/2𝜏 ௜ሻ𝑺෨௡∘ (17) 𝑯௡ାଵ௜≡𝑯෩௡௜൅e x p ሺെΔ𝑡/2𝜏 ௜ሻ𝑺෨௡ାଵ∘ (18) Using 𝑯෩௜, the algorithmic tangent operat or is obtained as follows. 𝑳௡ାଵ௔௟௚ൌ4 ቆ𝜕ଶ𝑈∘ 𝜕𝐽ଶ𝜕𝐽 𝜕𝑪⊗𝜕𝐽 𝜕𝑪൅𝜕𝑈∘ 𝜕𝐽𝜕ଶ𝐽 𝜕𝑪𝜕𝑪ቇ൅𝑔 ሺ∆𝑡/2 ሻ𝑳௡ାଵ∘തതതതതത ൅෍െ2 3𝑔௜𝐽௡ାଵିଶ/ଷቆDEV ௡ାଵ൫𝑯෩௡௜൯⊗𝑪௡ାଵିଵ൅𝑪௡ାଵିଵ ⊗D E V ௡ାଵ൫𝑯෩௡௜൯𝒏 𝒊ୀ𝟏 ൅൫ 𝑯෩௡௜:𝑪௡ାଵ൯൬𝑰𝑪೙శభషభ൅1 3𝑪௡ାଵିଵ⊗𝑪௡ାଵିଵ൰ቇ(19) A detailed derivation of the elastic tangent operator and algor ithmic tangent operator is presented in the literature [34]. 8 3. Adaptive affine homogenization in the absence of interfacial imperfection 3.1 Adaptive affine homogenization We propose the adaptive affine homogenization method which is d esigned to resolve a common drawback of existing methods, the inconsistency of the accumulated strain state and the concentration tensor. Our method utilizes the Euler forward algorithm in order to avoid an iterative process for numerical i ntegrations and to facilitate the convergence to an accurate solution with small time steps. We use the primary procedure of the incrementally affine method [25, 26] as a preliminary homogenization step to conside r inelastic stress change. In the process, each phase of the composite is given the same affi ne deformation to have the same stress. And then, macro deformation is assigned to remove the v irtually loaded affine deformation through the conventional homogenization process, wh ich resulted in the strain increments as below, ∆𝜺ଵൌ𝑨ଵఌ:ሼ∆𝛆ത൅𝑐 ଴𝕻:ሺ∆𝝈଴௥െ∆ 𝝈ଵ௥ሻሽ, ∆ 𝜺തൌ𝑐 ଴∆𝜺଴൅𝑐 ଵ∆𝜺ଵ (20) where 𝑐଴ and 𝑐ଵ are volume fraction of the matrix and particle, respectively, ∆𝜺ത, ∆𝜺଴, and ∆𝜺ଵ are the strain increments of the composite, matrix, and partic le at a given time step, respectively. In the present paper, subscripts 0 and 1 indicate state variables of the matrix and inclusion phase, respectively. Inelastic stress changes, i.e., stress relaxation of both phases due to viscosity, are shown as ∆𝝈଴௥ and ∆𝝈ଵ௥. 𝑨ଵఌ and 𝑩ଵఌ are strain concentra tion tensors which are defined in the linear elastic regime as follows, 𝑨ଵఌൌ𝑩ଵఌ:ሾ𝑐଴𝑰൅𝑐 ଵ𝑩ଵఌሿିଵwhere ∆𝜺 ଵൌ𝑨ଵఌ:∆ 𝜺ത (21) 𝑩ଵఌൌሾ𝑰൅𝕻 : ሺ𝑳ଵെ𝑳 ଴ሻሿିଵwhere ∆𝜺 ଵൌ𝑩ଵఌ:∆ 𝜺 ଴ (22) Here, 𝑰 and 𝕻 are the fourth-order identity tensor and Hill tensor [35]. 𝑳଴ and 𝑳ଵ a r e tangent operator of the matrix a nd particle, respectively. Deta iled explanation for the incrementally affine method can be found in previous papers [25 , 26]. In addition to the procedure prescribed by the incrementally af fine method, the adaptive affine method considers the modulus change in the defo rmed state. Because the modulus change over applied loading or time occurs in most nonl inear composites, our method 9 can be expanded to most types of LCC schemes of nonlinear compo site homogenization. We note that the modulus change dur ing each incremental deformatio n step leads to the mismatch between average strains over each phase and the strain concentr ation tensor. Moreover, such mismatch accumulates with the repeated incremental deformation procedures. Thus, in the adaptive affine method, a fictitious uniform strain ( ∆𝜺∗) is augmented to each phase of the composite in order to ensure equality, ∆𝜺ଵൌ𝑩ଵఌ:∆ 𝜺 ଴. After the augmentation of the fictitious uniform strain, the stress states of matrix and inclusions are written as below, 𝝈଴∗ൌ𝑳 ଴:∆ 𝜺∗൅𝝈 ଴ሺ𝑡௡ሻ (23) 𝝈ଵ∗ൌ𝑳 ଵ:∆ 𝜺∗൅𝝈 ଵሺ𝑡௡ሻ (24) The fictitious strain is then tuned to impose the consistency o f the strain tensor at the deformed state at 𝑡௡ାଵ, as below, 𝑳ଵିଵሺ𝑡௡ାଵሻ:𝝈ଵ∗ൌ𝑩ଵఌ:𝑳଴ିଵሺ𝑡௡ାଵሻ:𝝈଴∗. (25) Substituting Eq. (23) and Eq. (24) into Eq. (25), fictitious un iform strain ( ∆𝜺∗) is expressed as follows, ∆𝜺∗ൌሺ𝑰െ𝑩ଵఌሻିଵ:ൣ 𝑩ଵఌ:𝑳଴ିଵሺ𝑡௡ାଵሻ:𝝈଴ሺ𝑡௡ሻെ𝑳 ଵିଵሺ𝑡௡ାଵሻ:𝝈ଵሺ𝑡௡ሻ൧ (26) Afterwards, as in the incrementally affine method, an identical macro strain is applied to the composite in order to remove the fictitious deformation. Follow ing the procedure visualized in Fig. 1, the incremental strain cha nge of inclusion satisfying t he equality after the incremental deformation is expressed as follows, ∆𝛆ଵ௡∗ൌ𝑐 ଴𝑨ଵఌ:൫ 𝑳଴ିଵ:𝝈଴ሺ𝑡௡ሻെሼ𝑰൅𝕻 : ሺ𝑳ଵെ𝑳 ଴ሻሽ:𝑳ଵିଵ:𝝈ଵሺ𝑡௡ሻ൯. (27) 𝑳଴ and 𝑳ଵ are newly obtained tangent operators within [ 𝑡௡,𝑡௡ାଵ] which are different from the tangent operators within ሾ𝑡௡ିଵ,𝑡௡ሿ. We note that even though 𝜺଴ሺ𝑡௡ሻ and 𝜺ଵሺ𝑡௡ሻ a r e adjusted by ∆𝛆଴௡∗ and ∆𝛆ଵ௡∗ in the adaptive affine method, the macro strain of the composi te does not change during the process and satisfies the following condition, 𝑐଴∆𝜺଴௡∗൅𝑐 ଵ∆𝜺ଵ௡∗ൌ0. (28) 3.2 Application to visco-hyperelastic material 10 In the finite deformation scheme, there are several options of conjugate stress and strain variables and the corresponding ta ngent operators. As discussed in the literature [36], it is appropriate to use nominal stress rate, deformation gradient ra te, and nominal tangent operator because these stress and strain measures automatically satisfy the basic relation between macro measure and volume-weighted average scheme as follows, 𝑭ഥሶൌ𝑐 ଴𝑭ሶ଴൅𝑐 ଵ𝑭ሶଵ (29) 𝑷ഥሶ்ൌc ଴𝑷ሶ଴்൅c ଵ𝑷ሶଵ் (30) where 𝑷 ሶ்ൌ𝑫 :𝑭 ሶ, Here, 𝑷் refers to the nominal stress, and 𝑫 refers to the nominal tangent operator. The superposed dot represents the material time derivative. ሺ∙ሻതതത refers to volume averaged value. Our adaptive affine method is applied to study the mechanical r esponse of a highly nonlinear visco-hyperelastic material under finite deformation. The nonli nearity of visco-hyperelastic material has two primary origins. First, the tangent modulus ch anges with the strain as in a typical hyperelastic material. Second, the nonlinear viscoelast ic stress relaxation occurs as in a typical viscoelastic material. T he choice of a tangent operator is a critical issue to apply the adaptive affine method and the incrementally affine method to t he visco-hyperelastic composite. In section 2, the elastic (instantaneous) and algori thmic tangent operators were presented. The difference between elastic and algorithmic tange nt operators depends on time. The elastic tangent operator is defined at the time 𝑡௡ and thus the elastic tangent operator is identical to the instantaneous ta ngent operator. However, the a lgorithmic tangent operator is defined within the time-span ሾ𝑡௡,𝑡௡ାଵሿ. Since the stress relaxation occurs with deformation, the algorithmic tangent operator is usually softer than the ela stic tangent modulus mainly by the factor of 𝑔ሺΔt/2ሻ. From this point of view, it is reasonable to use the elastic tangent operator for the adaptive affine method while the algorithmic t angent operator for the incrementally affine method. Thus , in this paper, the following pseudo-code is devised. 1. Calculate elastic tangent ope rator of the matrix and inclusi on at 𝑡௡ 2. Convert all material tangent operators ( 𝑳଴∘ and 𝑳ଵ∘) to nominal tangent operators ( 𝑫଴∘ and 𝑫ଵ∘) as follows, 𝐷௄௜௅௝∘ൌ𝑆 ௄௅𝛿௜௝൅𝐹 ௜ெ𝐹௝ே𝐿௄ெ௅ே∘. (31) 11 3. Apply adaptive affine method to get adaptive deformation gra dients ( ∆𝑭଴௡∗ and ∆𝑭ଵ௡∗, Here, ∆𝜺௜௡∗ and 𝝈௜ in Eqs. (27, 28) are replaced to ∆𝑭௜௡∗ and 𝑃௜், respectively). 4. Calculate algorithmic tangent operators of the matrix and in clusion for time step [ 𝑡௡,𝑡௡ାଵ]. 5. Convert algorithmic material tangent operators ( 𝑳଴௔௟௚ and 𝑳ଵ௔௟௚) to algorithmic nominal tangent operators ( 𝑫଴௔௟௚ and 𝑫ଵ௔௟௚) 6. Apply incrementally affine method with macro deformation gra dient ( ∆𝑭ഥ). 7. Update deformation gradient increments ( ∆𝑭 ௜) and stress increments of each phase and go to step 1 for the next deformation increment. 3.3 Validation using FEM for perfect interfacial bonding case Our homogenization model is verif ied with the FEM calculation r esults, using a virtual visco- hyperelastic composite material in which finite deformation and viscoelastic relaxation affect its mechanical response significan tly. Commercial software ABAQ US and user-defined material (UMAT) are utilized for homogenization and FEM. We cho ose the Neo-Hookean model and prony series for the visco-hyperelastic matrix as fol lows, 𝜓ൌ𝐾 2ሺ𝐽െ1 ሻଶ൅𝐺 2൬𝐽ିଶ ଷtrሺ𝑪ሻെ3൰. (32) The inclusion also obeys a Neo-H ookean material, but it is much stiffer than the matrix and thus undergoes very small deformation. The material properties are summarized in Table 1. For the viscous response, only the deviatoric part is taken as a so urce of viscosity. We refer to values for normalized Prony series from a previous paper [17] a s shown in Table 2. The face centered cubic (FCC) structure i s adopted as the RVE model for FEM due to two reasons. First, the stress-strain curves of diff erent structures are not very d ifferent, and stress-strain curves of randomly distribution structure converge to that of FCC structu re with the numbe r of inclusions as shown in Fig. 2. Second, the FCC structure is symmetrical re ducing computational cost with good accuracy. Uniaxial cyclic loading tests with different str ain rates are conducted to compare our homogenization model and FEM results as shown in Fi g. 3. Bi-axial cyclic loading tests with different strain rates are conducted in Fig. 4. As shown in Fig. 3 and Fig. 4, 12 our homogenization model shows quite a good accuracy. Multi-axi al loading tests with a given loading profile are conducted for the different volume fraction of inclusions in Fig. 5. The results show that our model can be applied in a general loading case for visco-hyperelastic composite with good accuracy. 4. Adaptive affine homogenization with interfacial damage 4.1 Linear interfacial damage mod el-Effective inclusion method When interfacial damage exists between the matrix and inclusion , the overall mechanical response becomes softer than that of a perfect interfacial bond ing case. One of the widely adopted interfacial damage models is the interfacial spring mod el originally suggested by Qu et al [37, 38] in whic h the interfacial da mage is described by a linear spring la yer of vanishing thickness [39]. For the sake of mathematical simplicity, the mo del has been adopted in the mean-field homogenization scheme by either constructing the mod ified Eshelby tensor [39-41] or considering an effective inclus ion incorporating the interfa cial spring compliance [42, 43]. In this paper, we extend the effective inclusion method to fini te deformation regimes and implement the method in our homogenization process. Schematics of perfect bonding, interfacial spring model, and effective inclusion model are pre sented in Fig. 6. In contrast to the perfect bonding case , the interfacial spr ing model has a di splacement jump at the interface originated from the imperfect bonding or damage. The displaceme nt jump increases linearly with increasing interfacial tract ion because the virtual spring with compliance ( 𝛽) proportionally relates the displa cement jump to the interfacial traction. In this scheme, one has to consider several steps to calcu late the modified Eshelby ten sor and modified strain concentration tensor to implement the damage into the homogeniz ation process. Here, modified Eshelby tensor and modified strai n concentration tensor refer t o Eshelby tensor and strain concentration tensor in the presence of interfacial damage, res pectively. On the other hand, in the effective inclusion method, the displacement jump effect (i nterfacial damage) is fully considered by lowering the stiffness of the inclusion. After ob taining the damaged inclusion stiffness tensor, the user can regard it as a perfect bonding c ase. The detailed process to calculate the damaged inclusion stiffness tensor is shown below . The displacement jump is 13 expressed by the compliance of t he interfacial s pring as in the interfacial spring model. Δ𝑢௜ൌ𝛽 𝛥 𝑃 ௃௜𝑁௃ (33) where ∆𝒖 and β are the displacement jump and the compliance of the interfacia l spring, respectively. 𝑵 is the normal vector of the inte rface in the reference configu ration. Due to the displacement jump, macro defo rmation gradient increment ( Δ𝑭ഥ) is expressed as follows, Δ𝑭ഥൌ𝑐 ଴Δ𝑭 ଴൅𝑐 ଵΔ𝑭 ଵ൅𝑐 ଵ1 𝑉ଵනሺΔ𝒖 ⊗ 𝑵 ሻ𝑑𝑺 ப௏భ (34) ൌ𝑐 ଴𝛥𝑭 𝟎൅𝑐 ଵ𝛥𝑭 𝟏൅𝑐 ଵ𝛽 𝑎𝛥𝑷ଵ் (35) ൌ𝑐 ଴𝛥𝑭 ଴൅𝑐 ଵ൬𝑰൅𝛽 𝑎𝑫ଵ൰:𝑫ଵିଵ:𝛥 𝑷 ଵ (36) where 𝑎 represents the radius of the spherical inclusion. Compared to the perfect bonding case, the only difference is the third term in Eq. (34). As a result, the damaged inclusion stiffness tensor ( 𝑫ଵௗ) can be defined as follows, 𝑫ଵௗൌ𝑫 ଵ:൬𝑰൅𝛽 𝑎𝑫ଵ൰ିଵ (37) By replacing all 𝑫ଵ with 𝑫ଵௗ in the original homogenization process (Eq. (20) and Eq. (27)), the linear interfacial damage can be implemented on the homogen ization. 4.2 Validation using FEM for interfacial damage case The same material properties and structure as in Section 3.3 ar e used for the validation of interfacial damage cases. The interfacial damage is implemen ted by cohesive elements on the inclusion surface in FEM. Uniaxial tension tests are conduc ted with different values of 𝛽 as shown in Fig. 7. Our model shows a good match with FEM resul ts, but errors increase gradually with strain. We suspect that the errors originate fro m the difference in reference configuration between homogenizat ion and FEM. The virtual inter face spring is defined in the initial configuration in our mode l while it is redefined in eve ry current configuration in FEM. 5. Conclusion 14 In this study, we proposed the adaptive affine method which res olves a common limitation of existing homogenization methods, th e inconsistency between the strain concentration tensor and the accumulated strain states. The homogenization method sh ows a good match against FEM simulation results of RVEs under uniaxial or bi-axial cycli c loading conditions or even with multi-axial loading conditions with different loading prof iles. We extend our model to interfacial damage cases by adopting an effective inclusion met hod, which also shows a good match against FEM results for different values of linear interf acial compliance. Although this paper mainly focuses on visco-hyp erelastic composite, we believ e that the key idea of the adaptive affine scheme can be extended to various nonlinear hom ogenization models based on the LCC scheme where a dramatic change of tangent operators is involved. Acknowledgments The authors acknowledge support by the Basic Science Research P rogram (2019R1A2C4070690) and Creative Mate rials Discovery Program (20 16M3D1A1900038) through the National Resear ch Foundation of Korea (NRF). 15