Adaptive Affine Homogenization Method for Visco-hyperelastic Composites with Interfacial Damage

Despite intense research on the homogenization methods, it still is a challenging task to predict the nonlinear mechanical responses of visco-hyperelastic particulate-reinforced composites. In this work, we propose the adaptive affine method, a novel mean-field homogenization method designed to ensure the consistency of the accumulated strain state 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 method is also extended to predict the mechanical response in the presence of interfacial imperfections described by linear cohesive traction-separation laws. The analytic predictions are validated 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.


Introduction
Particle reinforced composites have been widely adopted to enhance various physical properties including the structural, thermal, thermoelectric, piezoelectric properties, as well as to gain multi-functionality [1][2][3][4]. For the industrial development of composites, it is a prerequisite to predict their properties in terms of the composition of their constituents and applied external mechanical, thermal, and electrical stimuli. In that regard, the homogenization theory is well established for the prediction of effective physical properties of composites under a linear response regime. However, there is still room for improvement when it comes to the nonlinear response originated from geometric or material nonlinearity. One typical example of particle reinforced composite in which both geometric and material nonlinearities play a critical role is elastomeric matrix composites. For example, solid propellant materials, which consist of rubbery matrix and oxidizer particles such as aluminum or ammonium perchlorate, have been extensively studied due to their importance in designing high precision rocket and missile propulsion systems [5][6][7][8][9][10]. Besides solid propellants, visco-hyperelastic elastomer composites are widely used in stretchable electronics, soft robotics, wearable systems, automobile industries, etc. Accordingly, a variety of homogenization frameworks that can appropriately account for such nonlinear responses have been suggested.
Many previous researchers have adopted mean-field homogenization (MFH) approaches to predict the effective properties of solid propellants [11][12][13][14][15][16]. For the linear elastic response, one can apply the well-established MFH approaches such 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 inclusion problem. However, for the prediction of inelastic behaviors, multiple methods were suggested to apply MFH and it remains controversial which frameworks provide numerically efficient and mathematically rigorous predictions. A frequency-domain homogenization method was established for linear viscoelastic composites in which the correspondence principle holds 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 procedure is followed for the MFH of elasto-viscoplastic composites with an integral affine formulation [18]. Besides being costly, the numerical inversion of the LCT by a collocation method might be inaccurate. Possible improvements have been studied by Rekik and Brenner [19].
On the other hand, significant efforts have been devoted to developing time-domain homogenization methods for the sake of mathematical simplicity and convenience. The linear comparison composite (LCC) concept was proposed to extend linear 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 formulations [20,21] or the linearization of the local constitutive models. Depending on the linearization scheme, secant or tangent operators are adopted [22][23][24]. The incrementally affine method was developed using the analogy between linearized elasto-viscoplasticity and thermoelasticity [25] and between linearized viscoelasticity-viscoplasticity and thermoelasticity [26]. However, the incrementally affine method tends to show stiff predictions in the viscoelastic regime. The prediction becomes stiffer in visco-hyperelastic material under finite deformation.
Alternatively, the incremental secant method was developed which combines the advantages of secant modulus and tangent modulus methods [27][28][29][30]. The method showed acceptable predictability for materials with elasto-plasticity or elasto-viscoplasticity, but its extension to viscoelastic-viscoplastic 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 : ; : ⊗

Visco-hyperelastic material modeling
In this paper, we adopt Simo's visco-hyperelastic model because of the following advantages [31][32][33]. First, the volumetric or deviatoric part can be selectively chosen as a source of viscosity because the formulation separates the volumetric 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 effects by using the time-temperature superposition principle.
Elastic stress is derived from strain energy density ( ∘ ) which is additively decomposed into volumetric part ( ∘ ) and deviatoric part ( ∘ ) as follows, Then, the second Piola-Kirchhoff stress ( ) is derived from strain energy density with a viscous response as follows, 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 description which is defined as follows, Viscoelastic properties are represented with the Prony series as follows, 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, The elastic tangent operator also can be considered as an instantaneous 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, is called an internal algorithmic variable which has an integral form. The detailed process to calculate is shown below.
is newly defined from to separate variables within the time step [ , ].
Using , the algorithmic tangent operator is obtained as follows.
A detailed derivation of the elastic tangent operator and algorithmic tangent operator is presented in the literature [34].

Adaptive affine homogenization
We propose the adaptive affine homogenization method which is designed 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 integrations 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 consider inelastic stress change. In the process, each phase of the composite is given the same affine deformation to have the same stress. And then, macro deformation is assigned to remove the virtually loaded affine deformation through the conventional homogenization process, which resulted in the strain increments as below, where and are volume fraction of the matrix and particle, respectively, ∆ , ∆ , and ∆ are the strain increments of the composite, matrix, and particle 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 concentration tensors which are defined in the linear elastic regime as follows, : where ∆ : ∆ Here, and are the fourth-order identity tensor and Hill tensor [35]. and are tangent operator of the matrix and particle, respectively. Detailed explanation for the incrementally affine method can be found in previous papers [25,26].
In addition to the procedure prescribed by the incrementally affine method, the The fictitious strain is then tuned to impose the consistency of the strain tensor at the deformed state at , as below, Afterwards, as in the incrementally affine method, an identical macro strain is applied to the composite in order to remove the fictitious deformation. Following the procedure visualized in and are newly obtained tangent operators within [ , ] which are different from the tangent operators within , . We note that even though and are adjusted by ∆ * and ∆ * in the adaptive affine method, the macro strain of the composite does not change during the process and satisfies the following condition,

Application to visco-hyperelastic material
In the finite deformation scheme, there are several options of conjugate stress and strain variables and the corresponding tangent operators. As discussed in the literature [36], it is appropriate to use nominal stress rate, deformation gradient rate, 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, 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 response of a highly nonlinear visco-hyperelastic material under finite deformation. The nonlinearity of visco-hyperelastic material has two primary origins. First, the tangent modulus changes with the strain as in a typical hyperelastic material. Second, the nonlinear viscoelastic stress relaxation occurs as in a typical viscoelastic material. The choice of a tangent operator is a critical issue to apply the adaptive affine method and the incrementally affine method to the visco-hyperelastic composite. In section 2, the elastic (instantaneous) and algorithmic tangent operators were presented. The difference between elastic and algorithmic tangent operators depends on time.
The elastic tangent operator is defined at the time and thus the elastic tangent operator is identical to the instantaneous tangent operator. However, the algorithmic tangent operator is defined within the time-span , . Since the stress relaxation occurs with deformation, the algorithmic tangent operator is usually softer than the elastic 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 tangent operator for the incrementally affine method. Thus, in this paper, the following pseudo-code is devised. 6. Apply incrementally affine method with macro deformation gradient (∆ ).

Update deformation gradient increments (∆ ) and stress increments of each phase and go
to step 1 for the next deformation increment.

Validation using FEM for perfect interfacial bonding case
Our homogenization model is verified with the FEM calculation results, using a virtual viscohyperelastic composite material in which finite deformation and viscoelastic relaxation affect its mechanical response significantly. Commercial software ABAQUS and user-defined material (UMAT) are utilized for homogenization and FEM. We choose the Neo-Hookean model and prony series for the visco-hyperelastic matrix as follows, The inclusion also obeys a Neo-Hookean 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 source of viscosity. We refer to values for normalized Prony series from a previous paper [17] as shown in Table 2. The face centered cubic (FCC) structure is adopted as the RVE model for FEM due to two reasons. First, the stress-strain curves of different structures are not very different, and stress-strain curves of randomly distribution structure converge to that of FCC structure with the number of inclusions as shown in Fig. 2. Second, the FCC structure is symmetrical reducing computational cost with good accuracy. Uniaxial cyclic loading tests with different strain rates are conducted to compare our homogenization model and FEM results as shown in Fig. 3. Bi-axial cyclic loading tests with different strain rates are conducted in Fig. 4. As shown in Fig. 3 and Fig. 4, our homogenization model shows quite a good accuracy. Multi-axial 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.

Linear interfacial damage model-Effective inclusion method
When interfacial damage exists between the matrix and inclusion, the overall mechanical response becomes softer than that of a perfect interfacial bonding case. One of the widely adopted interfacial damage models is the interfacial spring model originally suggested by Qu et al [37,38] in which the interfacial damage is described by a linear spring layer of vanishing thickness [39]. For the sake of mathematical simplicity, the model has been adopted in the mean-field homogenization scheme by either constructing the modified Eshelby tensor [39][40][41] or considering an effective inclusion incorporating the interfacial spring compliance [42,43].
In this paper, we extend the effective inclusion method to finite deformation regimes and implement the method in our homogenization process. Schematics of perfect bonding, interfacial spring model, and effective inclusion model are presented in Fig. 6 where ∆ and β are the displacement jump and the compliance of the interfacial spring, respectively. is the normal vector of the interface in the reference configuration. Due to the displacement jump, macro deformation gradient increment (Δ ) is expressed as follows, : : 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, : By replacing all with in the original homogenization process (Eq. (20) and Eq. (27)), the linear interfacial damage can be implemented on the homogenization.

Validation using FEM for interfacial damage case
The same material properties and structure as in Section 3.3 are used for the validation of interfacial damage cases. The interfacial damage is implemented by cohesive elements on the inclusion surface in FEM. Uniaxial tension tests are conducted with different values of as shown in Fig. 7. Our model shows a good match with FEM results, but errors increase gradually with strain. We suspect that the errors originate from the difference in reference configuration between homogenization and FEM. The virtual interface spring is defined in the initial configuration in our model while it is redefined in every current configuration in FEM.

Conclusion
In this study, we proposed the adaptive affine method which resolves a common limitation of existing homogenization methods, the inconsistency between the strain concentration tensor and the accumulated strain states. The homogenization method shows a good match against FEM simulation results of RVEs under uniaxial or bi-axial cyclic loading conditions or even with multi-axial loading conditions with different loading profiles. We extend our model to interfacial damage cases by adopting an effective inclusion method, which also shows a good match against FEM results for different values of linear interfacial compliance. Although this paper mainly focuses on visco-hyperelastic composite, we believe that the key idea of the adaptive affine scheme can be extended to various nonlinear homogenization models based on the LCC scheme where a dramatic change of tangent operators is involved.