The Effect of Evolving Damage on the Finite Strain Response of Inelastic and Viscoelastic Composites

A finite strain micromechanical model is generalized in order to incorporate the effect of evolving damage in the metallic and polymeric phases of unidirectional composites. As a result, it is possible to predict the response of composites with ductile and brittle phases undergoing large coupled inelastic-damage and viscoelastic-damage deformations, respectively. For inelastic composites, both finite strain elastoplastic (time-independent) and viscoplastic (time-dependent) behaviors are considered. The ductile phase exhibits initially a hyperelastic behavior which is followed by an inelastic one, and its analysis is based on the multiplicative split of its deformation gradient into elastic and inelastic parts. The embedded damage mechanisms and their evolutions are based on Gurson’s (which is suitable for the modeling of porous materials) and Lemaitre’s finite strain models. Similarly, the polymeric phase exhibits large viscoelastic deformations in which the damage evolves according to a suitable evolution law that depends on the amount of accumulated deformation. Evolving damage in hyperelastic materials can be analyzed as a special case by neglecting the viscous effects. The micromechanical analysis is based on the homogenization technique for periodic multiphase materials, which establishes the strong form of the Lagrangian equilibrium equations. These equations are implemented together with the interfacial and periodic boundary conditions, in conjunction with the current tangent tensor of the phase. As a result, the instantaneous strain concentration tensor that relates the local deformation gradient of the phase to the externally applied deformation gradient is established. This provides also the instantaneous effective stiffness tangent tensor of the composite as well as its current response. Results are given that exhibit the effect of damage on the initial yield surfaces, response and possible failure of the composite.

Abstract: A finite strain micromechanical model is generalized in order to incorporate the effect of evolving damage in the metallic and polymeric phases of unidirectional composites. As a result, it is possible to predict the response of composites with ductile and brittle phases undergoing large coupled inelastic-damage and viscoelastic-damage deformations, respectively. For inelastic composites, both finite strain elastoplastic (time-independent) and viscoplastic (time-dependent) behaviors are considered. The ductile phase exhibits initially a hyperelastic behavior which is followed by an inelastic one, and its analysis is based on the multiplicative split of its deformation gradient into elastic and inelastic parts. The embedded damage mechanisms and their evolutions are based on Gurson's (which is suitable for the modeling of porous materials) and Lemaitre's finite strain models. Similarly, the polymeric phase exhibits large viscoelastic deformations in which the damage evolves according to a suitable evolution law that depends on the amount of accumulated deformation. Evolving damage in hyperelastic materials can be analyzed as a special case by neglecting the viscous effects. The micromechanical analysis is based on the homogenization technique for periodic multiphase materials, which establishes the strong form of the Lagrangian equilibrium equations. These equations are implemented together with the interfacial and periodic boundary conditions, in conjunction with the current tangent tensor of the phase. As a result, the instantaneous strain concentration tensor that relates the local deformation gradient of the phase to the externally applied deformation gradient is established. This provides also the instantaneous effective stiffness tangent tensor of the composite as well as its current response. Results are given that exhibit the effect of damage on the initial yield surfaces, response and possible failure of the composite.
Keywords: large deformations; finite inelasticity; finite viscoelasticity; evolving damage; micromechanical analysis; finite strain high fidelity generalized method of cells

Introduction
In [1], a review of finite strain micromechanical analyses of multiphase materials have been presented. It was shown that it is possible to predict the microscopic (local) and macroscopic (global) response of composites undergoing large deformations in which the constituents in these composites can be modeled as hyperelastic, thermoelastic (based on entropic elasticity), viscoelastic (including quasilinear viscoelasticity (QLV) which is suitable for the modeling of biological tissues), thermoviscoelastic, rate-dependent thermoinelastic (viscoplastic) and rate-independent thermoinelastic (elastoplastic) materials. In all cases, the micromechanical analyses were based on the homogenization technique for periodic composites, and able to provide the composite's behavior in conjunction with the known properties of its constituents, their constitutive relations, detailed interaction and volume ratios. These analyses, referred to as High Fidelity Generalized Method of Cells (HFGMC), provide the instantaneous mechanical, thermal and inelastic concentration tensors that relate the local induced strain in the phase to the current externally applied strains and temperature. In addition, these micromechanical analyses yield the macroscopic constitutive equations of the multiphase composite in terms of its instantaneous stiffness (tangent) and thermal stress tensors. In anyone of these micromechanical analyses, the local field distribution among the various constituents of the composite can be also determined at any instant of loading.
Continuum damage considerations were not included in these investigation except for the prediction of the response of finite strain elastic composites where the Mullins effect can take place in the hyperelastic phase (Aboudi [2]). The purpose of the present paper is to include continuum damage mechanics considerations (see [3], [4], [5], [6], [7], for example) in finite strain elastoplastic, viscoplastic and viscoelastic materials any one of which can be a constituent in a multiphase composite. As a result, the presented generalized micromechanical analyses will enable the prediction of the behavior of composites with evolving damage in their phases.
In the case of elastoplasticity coupled with evolving damage, the model of Gurson [8] is employed as a phase in the composite. Gurson's model is capable to describe the behavior of homogenized porous materials in which the amount of porosity corresponds to the damage variable which evolves with the applied loading from a non-zero initial value and controlled by the plastic flow. In the present paper, we adopt the two-dimensional Gurson's model according to which the macroscopic response of a porous material is established by analyzing the field of a hollow cylinder that is subjected to a remote loading. Hence the response of the elastoplastic material in this case exhibits a transverse isotropy where the behavior in the axial direction along which the cylinder is oriented is different from the behavior in the transverse direction. The HFGMC model on the other hand can also predict the behavior of a finite strain porous elastoplastic material (which is just a special case of a composite material in which the stiffness of the second phase is negligibly small). Hence it is possible to compare the prediction of the the initial yield surfaces and response that are provided by the two independent models. Gurson's model is also employed to represent a phase in HFGMC method. This latter modeling implies the accommodation of inelastic anisotropic constituent in a composite (microscale level) to obtain its overall inelastic anisotropic behavior (macroscale level). The considered other finite strain elastoplasticity coupled with damage model in the present investigation, is the one presented by Lemaitre [9], [10]. Here the damage variable in the material evolves with applied loading from any initial value, and the evolution is controlled by the of plastic flow rule. Both Gurson's and Lemaitre's models were originally formulated in the framework of small strains, and their extension to finite strain was presented by de Souza Neto et al. [11]. In both models, the multiplicative decomposition of the deformation gradient is employed followed by the integration of the governing equations by means of the return mapping algorithm [12], [13], in conjunction with the logarithmic strain and the backward exponential approximation (Weber and Anand [14], Eterovic and Bathe [15], Cuitino and Ortiz [16], Simo [17]).
The analysis of finite strain viscoplasticity coupled with evolving damage in a material follows the same methodology of finite elastoplasticity. Here however the consistency parameter which appears in the flow rule is not one of the unknowns but it is prescribed in advance in terms of the other unknown field variables [12] (cf. Perzyna's [18]). Both the finite strain Gurson's and Lemaitre's models are presently extended and applied to investigate the behavior of composites that consist of viscoplastic coupled with evolving damage constituents.
A finite strain viscoelasticity theory was presented by Reese and Govindjee [19] which, in particular, allows large deviations away from the equilibrium state. This is in contrast to finite linear viscoelasticity theories, see Simo [20] and Holzapfel [21] for example and references cited in [19], which although account for large deformations, but restrict the formulation to states close to thermodynamic equilibrium by choosing linear evolution laws for the internal variables. In [19], the multiplicative split of the deformation gradient is also used, and the free energy of the viscoelastic material is expressed by the Ogden's representation (Ogden [22], Holzapfel [21]). The integration of the viscoelastic evolution equations is performed in conjunction with the logarithmic strain and the backward exponential approximation which yield a system of nonlinear equation. Extensive comparisons between the results obtained by the theories of Reese and Govindjee [19] and Simo [20] are given by Govindjee and Reese [23]. In the finite viscoelasticity theory of Reese and Govindjee [19] damage is not taken into account. In the present investigation, evolving damage in the finite viscoelastic material is included by adopting the derivation of Lin and Schomburg [24] and Miehe and Keck [25] according to which the rate of damage depends upon the kinematic arc-length. By neglecting the viscous effects, the special case of a hyperelastic material with evolving damage is obtained (it should me noted that Gurson's and Lemaitre's models do not reduce to such a special case because in the absence of inelasticity, damage does not evolve).
In all cases mentioned above the constitutive equations of the finite strain inelastic and viscoelastic monolithic materials are ultimately given in terms of the first tangent tensors that relate the rates of the first Piola-Kirchhoff stress and the deformation gradient tensors. The established first tangent tensors of the various constituents are subsequently employed in the finite strain HFGMC method for the prediction of the behavior of the resulting multiphase composites.
After briefly presenting the aforementioned derivations of inelastic and viscoelastic theories coupled with evolving damage together with the corresponding computational procedures, applications are given which show the behavior of the various types of composites in various circumstances. The paper is concluded by some suggestions for further generalizations in a future research.

Finite Strain Coupled Elastoplasticity-Damage Models of Monolithic Materials
Let X and x denote the location of a point in the material with respect to the initial (Lagrangian) and current systems of coordinates, respectively, and t is the time. In terms of the local deformation gradient tensor F(X, t), dx = F(X, t)dX. The finite plasticity theory that is employed in the present investigation is based on the introduction of a stress-free intermediate configuration and a multiplicative decomposition of the local deformation gradient F(X, t) in the form where F p (X, t) and F e (X, t) are the deformation gradient tensors from the initial to the intermediate and from the intermediate to the current configuration, respectively. The corresponding right Cauchy-Green tensors are given by where superscript T denotes the transpose operation. The left Cauchy-Green tensors B and B e are defined by The logarithmic elastic strain e is defined by Two types of energy functions that model the monolithic phase in the composite are considered in the present investigation. The first one is based on the macroscopic yield function of homogenized porous materials that was established by Gurson [8]. Here, the damage variable corresponds to the amount of porosity and its evolution describes its growth with the applied loading. The second energy function is due to Lemaitre [9], [10] which includes an evolving damage variable. Both models were originally formulated in the framework of infinitesimal deformations, and their extension to finite strain analysis was presented by de Souza Neto et al. [11]. These models are briefly described in the following.

Gurson's finite strain elastoplastic porous material model
The Hencky isotropic strain energy function per unit reference volume is defined by where the fourth-order tensor h given by where λ and µ are material parameters, and I and I 4 being the identity second and fourth-order tensors, respectively. The Kirchhoff stress τ is given by The macroscopic yield function for homogenized porous materials that was established by Gurson [8] in which the long cylindrical pores (voids) are oriented in the 1-direction is given by where dev [τ ] is the deviator of Kirchhoff stress tensor, D denotes the porosity volume fraction (which represents the amount of damage) and σ y is the function that describes the hardening law of the material. For isotropic hardening it is given by where Y 0 is the yield stress in simple tension and K(R) describes the isotropic hardening law. For linear hardening: K(R) = HR. The plastic flow rule is determined from Φ by employing the relation [12]: where L v is the Lie derivative of B e and γ is the consistency parameter. Here and in the followingQ means the derivative of a quantity Q with respect to time or quasi time t. The Lie derivative of B e can be expressed in the form [13] L It is worth mentioning that the flow rule (10) is identical to the one used in the book by Bonet and Wood [26], and the von Mises plasticity is readily obtained from Equation 8 by setting D = 0. The finite strain plastic flow rule of Simo and Hughes [13] that was employed by Aboudi [27] however is given, apart from a numerical coefficient, by Equation 10 but with [B e ] −1 on the left-hand-side replaced by (trace [B e ]) −1 . It is worth mentioning that both the present and Simo and Hughes [13] flow rules provide almost identical responses of the ductile material that is specified in the following when it is subjected to uniaxial stress loading. The responses significantly differ however for uniaxial strain loading although yielding occurs at the same strain. The rate of hardening is determined from the relation [11] whereas the evolution law for damage (porosity) is given bẏ Equations 10-13 together with the condition that Φ = 0 are solved by employing the return mapping algorithm [11]. It should be emphasized that by means of the backward exponential approximation (Weber and Anand [14], Eterovic and Bathe [15], Cuitino and Ortiz [16], Simo [17]) the plastic flow rule that is given by Equation 10 can be reduced to the simple update expression e n+1 = e trial n+1 − ∆γ ∂Φ ∂τ n+1 (14) where e trial n+1 at the present n + 1 increment is computed according to Equation 4 while employing B e n of the previous increment. The fourth-order first tangent tensor R that is needed in the micromechanical analysis is defined bẏ where T is the first Piola-Kirchhoff stress tensor. It is related to the Kirchhoff stress tensor τ according to: T = F −1 τ . The first tangent tensor is determined according to : : : where  (17) and δ ij is the Kronecker delta. The tensor ∂τ /∂ e trial n+1 is determined by the differentiation of the set of equations that appear in the return mapping algorithm, see de Souza Neto et al. [11] for details.
Remark: It should be emphasized that a critical check to the validity of the established first tangent tensor R (here and in all the following cases) is that the values of the first Piola-Kirchhoff T that are obtained by the integration of Equation 15 must be identical to the values of T = F −1 τ that are computed directly from Equation 7 in which R is not involved.

Lemaitre finite strain elastoplastic damage model
The extension of Lemaitre elastoplastic model [9], [10] that includes an evolving damage to large deformations was presented by de Souza Neto et al. [11]. It is based on Hencky isotropic elastic-damage free energy function per unit reference volume (cf. Equation 5): where as before e denotes the logarithmic elastic strain and D represents the damage variable such that 0 ≤ D ≤ 1. The Kirchhoff stress is determined according to, cf. Equation 7, The thermodynamical force Y which is conjugate to the damage variable D is given by By assuming isotropic hardening, the yield function Φ of this model is given by where J 2 (τ ) = Finally, the dissipation function is taken in the form where r and s are material constants. This function provides the evolution law of damage in the forṁ Equations 19-24 together with the condition that Φ = 0 are solved by employing the return mapping algorithm [11]. The backward exponential approximation reduces the flow rule which is given by Equation 10 to the simple relation (14). The final constitutive relation of Lemaitre's finite strain elastoplasticdamage model can be expressed in the form shown by Equation 15 where the fourth-order first tangent tensor R is determined by Equation 16. As in the Gurson's model, the tensor ∂τ /∂ e trial n+1 that appears in Equation 16 is determined by the differentiation of the set of equations that appear in the return mapping algorithm.
Finally, it is worth mentioning that unlike Gurson's model, the value of damage D evolves, in the framework of Lemaitre's model, from D = 0 according to Equation 24. In Gurson's model however, a damage variable whose initial value is equal to zero does not evolve, see Equation 13.

Finite Strain Coupled Viscoplasticity-Damage Models of Monolithic Materials
The previously discussed Gurson and Lemaitre damage models were based on finite strain elastoplasticity (time-independent). Both models can be generalized to exhibit viscoplastic (time-dependent) behavior. In this situation the time derivative of the consistency parameterγ that appears in Equation 10 is not one of the unknowns that need to be determined as in the elastoplasticity case. Following the methodology that was presented in [12],γ in the coupled viscoplastic-damage Gurson's model can be established in the formγ where the yield function Φ is given by Equation , η and are parameters that characterize the viscoplastic material, and H(Φ) is the Heaviside unit step function.
For Lemaitre coupled viscoplastic-damage model,γ is given bẏ where Φ is given by Equation 21. In the limit when η → 0 and → 0, the viscoplastic model recovers the perfectly plastic von Mises model [12]. It should be noted thatγ in Perzyna's [18] viscoplasticity model has a slightly different form from the above representation. For the Lemaitre model, for example, Equation 26 takes the formγ The final constitutive law of both Gurson and Lemaitre finite strain coupled viscoplasticity-damage models has the form given by Equation 15 where the fourth-order first tangent tensor R is given by Equation 16. In this latter equation, the instantaneous fourth-order tensor ∂τ /∂ e trial n+1 is determined as follows. From Equation 14 we obtain that It follows from Equation 7 that for Gurson's viscoplastic model where the following definitions are used Let us multiply h : ∆ vp in Equation 28 by (cf. Boyd and Allen [28]) This readily yields that It follows that the requested instantaneous tangent tensor for the viscoplastic Gurson's model is of the form which can be employed in Equation 16. For Lemaitre's viscoplastic model, Equation 19 should be employed. Hence Let us denote Hence the instantaneous tangent tensor for this model is given by which can be employed in Equation 16.

Finite Strain Coupled Viscoelasticity-Damage Model of Monolithic Materials
In the present section we briefly present the constitutive behavior of finite strain viscoelastic materials that exhibit evolving damage. The presentation follows the papers of Reese and Govindjee [19] where no damage is accounted to and Lin and Schomburg [24] where evolving damage is included. It should be mentioned that the present viscoelastic modeling allows large deviations from the thermodynamic equilibrium state and therefore is referred to as finite viscoelasticity (in contrast to finite linear viscoelasticity where small deviations from the equilibrium state is assumed).
The multiplicative decomposition of the deformation gradient is given this time by where F v is the counterpart of F p in Equation 1. The modeling that is presented herein is based on a single Maxwell and elastic elements, but it can be extended to include several Maxwellian elements. The total free energy per unit reference volume is decomposed into equilibrium (EQ) which represents the strain energy of the elastic element and a nonequilibrium (NEQ) part that accounts for the Maxwell element: where ψ EQ 0 and ψ N EQ 0 are referred to as the effective free energy of the undamaged material, and D denotes the amount of damage such that 0 ≤ D ≤ 1.
The resulting Kirchhoff stresses are given by where C e = F eT F e and τ EQ 0 , τ N EQ 0 correspond to the effective Kirchhoff stresses of the undamaged material.
Let the left Cauchy-Green tensor B be represented in terms of its eigenvalues: With J = det F , the volume preserving tensorB = J −2/3 B can be accordingly represented in the form The finite strain elastic contribution can be modeled by the Ogden's compressible material representation (Ogden [22], Holzapfel [21]) as follows where K e is the elastic bulk modulus and and µ e r and α e r are material parameters of the elastic element. For Maxwell's element, the free energy is represented by [19]: where and The evolution equation in the present finite viscoelastic case that replaces Equation 10 is given by where η D and η V are the deviatoric and volumetric viscosities, respectively, and presently (cf. Equation  11) By employing the exponential mapping algorithm, Equation 47 is reduced to, cf. Equation 14, where the principal values of the logarithmic strain e A are given by, cf. Equation 4, e A = 1/2 log(b e A ) and ∆t is the time increment between the current and previous step. Equation 49 forms a system of coupled nonlinear equations in the three unknowns: e A , A = 1, 2, 3. It readily provides, cf. Equation 27, that The rate of damage evolution is given by Lin and Schomburg [24] and Miehe and Keck [25] in the formḊ with the saturation value and In these relations, η dam , D ∞ 0 and α dam are material parameters. The fourth-order first tangent tensor R N EQ is determined by employing Equation 16. In this equation the expression ∂τ N EQ /∂ e trial n+1 is determined as follows. From Equation 39 the following expression can be established Let the second-order tensor M be defined by With the explicit components of M are given by [19]: In conjunction with Equation 50, we obtain that Let ∆ e A and ∆ ved A be defined by Therefore Equation 62 can be represented by By following the previous derivation for the establishment of the tangent tensor in the viscoplastic case, the following expression can be obtained With τ N EQ A , A = 1, 2, 3, and their derivatives given by Equations 39, 46 and 65, the derivative fourthorder tensor function ∂τ N EQ /∂ e trial n+1 can be established by employing the explicit expressions given by [11]. Similar treatment applies to τ EQ which leads to the establishment of ∂τ EQ /∂ e n+1 and R EQ . The first tangent tensor of the viscoelastic material is given by R = R N EQ + R EQ . The special case of a hyperelastic material coupled with evolving damage can be obtained by disregarding the viscous effects so that ψ N EQ = 0.

Finite Strain Micromechanical Analysis
Finite strain HFGMC micromechanical analyses for the establishment of the macroscopic constitutive equations of various types of composites with doubly periodic microstructure have been previously reviewed by Aboudi [1]. These micromechanical analyses are based on the homogenization technique in which a repeating unit cell of the periodic composite can identified. This repeating unit cell represents the periodic composite in which the double periodicity is taken in the transverse 2 − 3 plane, so that the axial 1-direction corresponds to the continuous direction (for a fiber-reinforced material, for example, the 1-direction coincides with the fibers orientation). In the framework of these HFGMC micromechanical models, the displacements are asymptotically expanded and the repeating unit cell is discretized. This is followed by imposing the equilibrium equations, the displacement and traction interfacial conditions as well as the conditions that ensure that the displacements and tractions are periodic across the repeating unit cell. In particular, the imposition of the equilibrium equations provide the strong form of the Lagrangian equilibrium conditions of the homogenization theory that must be satisfied. The resulting homogenization derivation establishes the deformation concentration tensor A(Y), where Y are the local Lagrangian system of coordinates with respect to which field variables in the repeating unit cell are characterized. This tensor relates the rate of the local deformation gradient gradientḞ(Y) at a material point Y within the repeating unit cell to the externally applied deformation gradient rateḞ in the form: It follows from Equation 66 and in conjunction with (15) that the local stress rate at this point is given byṪ Hence the resulting macroscopic constitutive equation for the multiphase composite undergoing large deformation is given byṪ where R * is the instantaneous effective stiffness (first tangent) tensor of the multiphase composite. It can be expressed in terms of the first tangent tensors of the constituents R(Y) and the established deformation concentration tensor A(Y) in the form where V Y is the volume of the repeating unit cell. More details can be found in the aforementioned reference [1]. It should be noted that the current value of R * of the composite is affected by the current value of damage variable through the instantaneous value of tangent tensors R(Y) of the finite strain constituents.
The finite strain HFGMC micromechanical model predictions were assessed verified by comparison with analytical and numerical large deformation solutions by Aboudi and Pindera [29], Aboudi [27] and Aboudi [2] for elastic, elastoplastic and elastic composite in which the Mullins damage effect is incorporated with the hyperelastic constituents, respectively.

Computational Procedures
In the following, the computational procedures for the determination of the finite strain response of the inelastic and viscoelastic composites are described.

Finite strain inelastic composites
At time t n , the local deformation gradient tensor F n and the elastic left Cauchy-Green tensor B e n have been already established in the inelastic phase of the composite. In addition, F n+1 is assumed to be known at time t n+1 = t n + ∆t.
(1) Compute the local trial logarithmic strain e trial n+1 as follows: (2) With e trial n+1 , compute τ trial n+1 from Equation 7 and solve the nonlinear system of equations that consists of Equation 14, the appropriate yield condition of the elastoplastic constituent, the evolution equations of the internal and damage variables in the constituent. In the case of a composite with a viscoplastic phase, the yield conditions of the constituent are not part of the system and, as previously mentioned, γ is not one of the unknowns. The solution of this system of nonlinear equations provides the current variables in the inelastic constituent of the composite at time step n + 1.
(3) The instantaneous first tangent tensor R of the constituent can be determined from Equation 16. (4) With the established local value of R in the inelastic phase, it is possible to proceed with the micromechanical analysis and compute the current concentration tensor A that appears in Equation 66, see [1]. Consequently, Equation 69 can be employed to determine the instantaneous effective tangent tensor R * of the composite. (5) The rate of the externally applied deformation gradient˙F is determined in accordance with the prescribed type of loading. For a uniaxial deformation in the 1-direction, for example,˙F 11 is known while all other components of˙F are zero. Hence, it is possible to compute the stress rate˙T using Equation 68. In addition, by integrating the resulting global stress rate,T is obtained. (6) Once˙F has been determined, it is possible to compute the rates of the local deformation gradients from Equation 66. The integration of the latter would yield the local deformation gradients in the constituent to be used in the next time step. (7) If, on the other hand, a uniaxial stress loading is applied on the composite such that˙F 11 only is known together withT the components of which are zero in the other directions, an iterative procedure is needed to determine the other components of˙F from these conditions.

Finite strain viscoelastic composites
At time t n , the local deformation gradient F n and the left Cauchy-Green deformation tensors B e n have been already established, and F n+1 is assumed to be known.
(1) From F n and F n−1 , the local right Cauchy-Green deformation tensors C n and C n−1 can be readily determined. Henceż in Equation 52 can be determined, from whichḊ that is given by Equation 51 can be integrated to provide the damage variable at time t n+1 : Equation 54 yields , see [1]. Consequently, Equation 69 can be employed to determine the tangent tensor R * . (7) The rate of the externally applied deformation gradient˙F is determined in accordance with the prescribed type of loading. For a uniaxial deformation in the 1-direction, for example,˙F 11 is known while all other components of˙F are zero. Hence, it is possible to compute the stress rate˙T using Equation 68. In addition, by integrating the resulting global stress rate,T is obtained. (8) Once˙F has been determined, it is possible to compute the rates of the local deformation gradients from Equation 66. The integration of the latter would yield the local deformation gradients to be used in the next time step. (9) If, on the other hand, a uniaxial stress loading is applied on the composite such that˙F 11 only is known together withT the components of which are zero in the other directions, an iterative procedure is needed to determine the other components of˙F from these conditions.

Applications
In this section, applications are given for the various models with evolving damage. For the inelastic constituent we consider a ductile material with linear hardening whose properties are given in Table 1 (In the range of small deformations these parameters correspond to the characterization of the aluminum alloy 2024-T4). The parameters of the finite viscoelastic material, Equation 44, are given in Table 2 together with η D and 1/η V = 0 (assuming elastic bulk deformations). The parameters E, ν, Y 0 and H denote the Young's modulus, Poisson's ratio, yield stress in simple tension and linear hardening slope. Table 2. Material parameters of the viscoelastic material (Reese and Govindjee [19]).
The parameters µ v p and α v p , p = 1, 2, 2 are the Ogden's material constants, K v is the bulk modulus and η D is the viscoelastic constant of the viscous element with 1/η v = 0.

Application of the finite strain Gurson's coupled elastoplastic-damage model
Gurson's coupled elastoplastic-damage two-dimensional model that was previously presented corresponds to the modeling of the effective behavior of a porous material in which the long cylindrical pores are oriented in the 1-direction. Consequently, it should be interesting to compare the predictions of the Gurson's model with the corresponding ones obtained from HFGMC micromechanical analysis. To this end, the doubly periodic HFGMC in which one of the phases has zero stiffness is considered. the volume fraction of this phase is taken to be equal to the porosity D in Equation 8. We start our investigation by comparing the initial yield surfaces as predicted by the Gurson's and HFGMC models.
For the Gurson's model the initial yielding is determined from Equation 8 which provides The von Mises yielding criterion is obtained from this equation by setting D = 0. Suppose that in the Gurson's model the initial yield surfaceT 22 againstT 33 is requested, whereT = F −1 τ being the first Piola-Kirchhoff stress tensor. Here all componentsT ij are equal to zero exceptT 22 andT 33 that may take any combination. Hence let us represent these components in the form where A is the radial distance of a point located on the initial yield surface envelopeT 22 -T 33 and θ is the corresponding polar angle. The substitution of these values in Equation 73, shows that for a given angle θ, the radial distance of A/Y 0 of this point is given by root of the transcendental equation where in Equation 73, τ 22 = F 22T22 and τ 33 = F 33T33 have been substituted. Since in typical metallic materials yielding takes place at very small strains (in aluminum, for example, yielding in simple tension accrues at a strain of about 0.4 percent), it can be practically assumed that F 22 = F 33 ≈ 1. Other types of initial yield surfaces can be generated in the same manner. The initial yield surfaces that are predicted by the HFGMC model, on the other hand, can be generated by establishing the instantaneous stress concentration tensor B(Y ) which relates the rate of the local Piola-Kirchhoff stress tensorṪ (Y ) to the externally applied stress rate˙T , namelẏ From Equations 66-68 this tensor can be shown to be given by Corresponding to the above discussion, suppose that the initial yield envelopeT 22  It can be clearly observed that a fair correspondence between the two models exist. Next, Figure 2 shows the initial yield surfaces in theT 11 -T 22 plane. Since the voids in the Gurson and HFGMC models are oriented in the 1-direction, the symmetry which can be observed in Figure 1 in the 2 and 3-directions does not exist anymore in the 1 and 2-direction in Figure 2. Here too, fair agreement between the two models can be observed. Finally, let us examine the initial yield envelopes for the loadingT 11 andT 22 =T 33 . The resulting envelopes are shown in Figure 3 for the above three values of porosity. It should be noted that the simple von Mises criterion (for yielding of homogeneous ductile materials) reduces in this case to two parallel straight linesT 11 −T 22 = ±Y 0 the first one of which passes the points: (T 11 /Y 0 = 1,T 22 /Y 0 = 0) and (T 11 /Y 0 = 0,T 22 /Y 0 = −1), whereas the second one passes through the points: (T 11 /Y 0 = −1,T 22 /Y 0 = 0) and (T 11 /Y 0 = 0,T 22 /Y 0 = 1), with the expected result that yielding will not occur for stress values between these two lines. Here too, the correspondence between the two models is quite good. In conclusion, these three figures show that the simple Gurson's model is quite reliable for the prediction of initial yield surfaces of porous materials.
Thus far the initial yield envelopes that are predicted by the two models have been investigated. The next investigation is concerned with the responses of porous materials that are obtained by the Gurson and HFGMC model. To this end, let us consider the ductile material that is characterized by Table 1. The response of this material with three values of porosity: D = 0.05, 0.25 and 0.4, is shown in Figure 4. This figure shows the average uniaxial stress responsesT 11 −F 11 to loading in the 1-direction (i.e., in the direction of which the pores are oriented) as predicted by the Gurson and HFGMC model. Also shown for reference is the uniaxial stress response of the bulk material (i.e., with zero porosity D = 0). The graphs show that good agreement between the predictions of Gurson and HFGMC model exists. The effect of porosity (damage) is clearly exhibited by comparison with the homogeneous material behavior.
The responsesT 22 −F 22 as predicted by the two models to a uniaxial stress loading of the porous material in the transverse 2-direction are shown in Figure 5 for three values of porosity (the response of the homogeneous material is included for comparison). Here too, very good agreement between the Gurson  and HFGMC model can be clearly observed. It should be noted that a careful comparison between the flow stress levels of Figures 4 and 5 reveals that the axial 1-direction along which the pore extends is relatively stronger (exhibiting higher stresses) than the transverse 2-direction. This observation is consistent with Figure 2 which shows that the porous material yields earlier when loaded in the transverse direction as compared to a loading in the axial direction. In Figures 4 and 5, the porous material behavior was shown in various cases at every one of which the prescribed value of porosity D was held constant. In the framework of Gurson's model the porosity (damage) can evolve as the plasticity develops, see Equation 13. Hence, it should be interesting to allow the damage, in the framework of Gurson's model, to evolve from a certain initial value D i to a final value D f which is determined by the model when the loading is terminated. The resulting response can be compared with HFGMC prediction in which the volume fraction of the pores is equal to D i and D f .  As a final investigation of the effect of evolving damage in the Gurson's model, we consider again, in the framework of the HFGMC, a composite material whose matrix is characterized by Table 1 while the second phase is a pore. Let D, as before, denotes the amount of porosity of this composite. As the loading of the composite (porous) material increases, damage evolves in the ductile matrix the amount of which is denoted by D m . Figure 7 shows the uniaxial stress responseT 22

Application of the finite strain Lemaitre's coupled elastoplastic-damage model
In the present subsection, a unidirectional composite that consists of a rubber-like hyperelastic matrix reinforced by ductile fibers whose properties are given in Table 1 is considered. The metallic fibers are modeled by Lemaitre's coupled elastoplastic-damage representation that was discussed in subsection 2.2, where the values of the parameters r and s are: r = 4.5M P a, s = 1. Figure 8(a) shows the uniaxial stress response response of the monolithic ductile material in the absence of any damage (i.e., D ≡ 0) and the corresponding case in which the damage evolves from D = 0 according to Equation 24. The response in the latter case is terminated when the damage approaches its final value D = 1. The evolution of the damage is exhibited in Figure 8(b) which shows its gradual increase with the applied loading. . Comparisons between the uniaxial stress response to loading in the axial 1direction of the porous ductile elastoplastic material whose matrix properties are given in Table 1   The rubber-like matrix is modeled as a hyperelastic compressible neo-Hookean material whose strain energy function is given by (Bonet and Wood [26]) where I 1 is the first invariant of the right Cauchy-Green deformation tensor C, J = det F and λ and µ are material constants. The second Piola-kirchhoff stress tensor S is obtained from and the constitutive relations of this material can be established in the form given by Equation 15, where the instantaneous first tangent tensor R of the material is given by  Comparisons between the uniaxial stress response to loading in the transverse 2-direction of the porous ductile elastoplastic material whose matrix properties are given in Table 1  Consider this metal/rubber-like composite in which the continuous metallic fibers are oriented in the axial 1-direction. The volume fraction of the fibers is v f = 0.1. Let the composite be subjected to an off-axis uniaxial stress loading. Here, the fibers which are oriented in the 1-direction, are rotated around the 3-direction by an angle φ. As a result, a new system of coordinates (X, Y, Z) is obtained such that Z = X 3 . The uniaxial stress loading is applied in the X-direction which is at angle φ with respect to the fibers direction. Referring to this new system of coordinates, the composite is loaded by the application of the deformation gradient F XX , and all components of the first Piola-Kirchhoff stress tensor T X , referred to the new coordinate system, are equal to zero except T XX . In particular, φ = 0 • and 90 • correspond to longitudinal and transverse uniaxial stress loading, respectively.
The locations of the initial yielding of this composite can be determined by generating its uniaxial stress response at various off-axis angles φ and detecting the stress at which yielding of the metallic fibers takes place. Figure 9 shows the off-axis stress-displacement gradient response of the composite in the absence of damage in the fiber phase (D ≡ 0). The locations of yielding points are indicated by the arrows. It should be noted that the earliest yielding is obtained when the composite is loaded in the axial direction (φ = 0 • ), whereas loading in the transverse direction (φ = 90 • ) caused yielding at a later Figure 6. Comparisons between the uniaxial stress response to loading in the transverse 2-direction of the porous ductile elastoplastic material whose matrix properties are given in Table 1   stage. Among the six off-axis angles at which the response of the composite is generated in the range of 1 ≤ F XX ≤ 3.5, the highest yield stress is obtained at φ = 60 • . Figure 10 presents that the metallic/rubber-like composite's response to uniaxial stress loading at various off-axis angles. Here, both cases in which the damage in the fiber phase is absent D ≡ 0 as well as when it is evolving from zero are shown. In the latter case, the dependence of damage evolutions on the applied deformation gradient are also shown. For loading in the axial direction φ = 0 • , the damage variable increases with applied loading until it reaches its final value of D = 1. For the off-axis angle φ = 20 • , on the other hand, the amount of evolving damage in the fiber phase is very small, as a result of which the difference between the responses in the absence and presence of damage is indistinguishable. It should be noted that the for φ = 0 • , 10 • and 15 • , the computations were terminated when the damage variables reached a certain value after which it jumped in the next step of computations to D = 1.  It is possible of course to increase the amplitude of applied loading in Figure 10 above the value of F XX = 1.5 to get appreciable values of damage in the fiber constituent for the off-axis loading case φ = 20 • . It is interesting however to observe the response and damage evolution when the metallic/rubberlike composite is subjected to a cyclic loading such that 0.5 ≤ F XX ≤ 1.5. Figure 11 presents the response of the composite and damage evolution in the fiber phase for a uniaxial stress cyclic loading for an off-axis angle of φ = 25 • . The figure clearly shows that the macroscopic stress T XX exhibits a cyclic behavior but the damage in the fiber constituent continue to increase while during unloading it retains its value. The figure indicates that the computations were terminated as the damage variable approached 1 after two cycles. Figure 8. (a) The uniaxial stress response of the monolithic metallic material whose properties are given by Table 1. Its constitutive relations are given by the Lemaitre elastoplastic model. The response is shown in the absence (D = 0) and evolving damage and, (b) the form of the damage evolution with the applied displacement gradient. (c) The uniaxial stress response of the monolithic hyperelastic matrix whose strain energy is given by Equation 79.

Application of the finite strain Gurson's coupled viscoplastic-damage model
Let us consider a porous viscoplastic material in which the properties of its matrix are given in Table 1. This porous material is modeled by Gurson's coupled viscoplastic-damage relations in conjunction with Equation 25 in which the viscosity and rate sensitivity parameters are: η = 1s and = 1, respectively. The applied rate of loading is 1s −1 . As in the elastoplastic case, the resulting response of the considered viscoplastic porous material can be compared with the prediction obtained from HFGMC in which the bulk matrix is described by Table 1 in conjunction with the above parameters and rate of loading. Figure 12 exhibits the response of the viscoplastic porous material with various amount of porosities as predicted by the Gurson's and HFGMC models. The porous material is uniaxially stress loaded in the axial 1-direction. Also shown is the corresponding behavior of the homogeneous viscoplastic material. In all cases there is no damage evolution. This figure is the viscoplastic counterpart of   Figure 13 where the porous material is uniaxially stress loaded in the transverse direction. Note that the scale of the plots in Figure 13 is different from that of Figure 12 or of the elastoplastic counterpart that is shown by Figure 5. The largest difference between Gurson's and HFGMC prediction is obtained for the lowest value of porosity (D = 0.05) and it decreases rapidly with increasing porosity.

Application of the finite strain Lemaitre's coupled viscoplastic-damage model
Consider the ductile material whose properties are given by Table 1. This material is represented by the Lemaitre's coupled viscoplastic-damage model, see Equation 26, in which the viscosity and rate sensitivity parameters are: η = 1s and = 1, respectively. In all cases shown herein the loading was applied at a rate of 1s −1 , and r = 4.5M P a, s = 1. Figure 14 shows a comparison between the response of the homogeneous ductile material when it is represented by Lemaitre's viscoplastic and elastoplastic coupled with damage models. The computations were terminated when the damage D approaches unity. This figure shows also the special case in which no damage takes place (i.e., D ≡ 0). As it is expected, the elastoplastic case exhibits lower stress values as compared to the viscoplastic behavior, but the evolution of damage is more rapid in the viscoplastic case.  Figure 12. Comparisons between the uniaxial stress response to loading in the axial 1direction of the porous viscoplastic material whose matrix properties are given in Table 1   The behavior of a unidirectional metal/rubber-like composite (with fiber volume fraction v f = 0.1) is shown in Figure 15 in which the metallic fibers (whose properties are given by Table 1) are represented by Lemaitre's viscoplastic and elastoplastic coupled with damage models, whereas the rubber-like matrix behavior is given by Equation 79. The responses shown in this figure are caused by the application of Figure 13. Comparisons between the uniaxial stress response to loading in the transverse 2-direction of the porous ductile viscoplastic material whose matrix properties are given in Table 1 Figure 14. A comparison between the response and damage evolution of the monolithic ductile material, whose properties were specified by Table 1, when it is modeled by Lemaitre viscoplastic (VP) and elastoplastic (EP) equations. Also shown are the corresponding special cases when the damage is ignored (D ≡ 0). (a) Stress-deformation response and, (b) damage evolution. uniaxial stress loadings F XX applied at off-axis angles φ = 0 • and 10 • with respect to the fibers. The figure shows also the damage evolution in the metallic phase and the response in the special case in which the damage in the fibers is ignored. For the axial loading case (φ = 0 • ) the elastoplastic and viscoplastic response and damage evolution are quite similar. The graphs show however that in contradistinction to the elastoplastic case (Figure 10c,d), damage evolution in the viscoplastic case when φ = 10 • is relatively small as compared to the elastoplastic prediction in which the damage approaches unity. As a result, the stress-deformation curve in the viscoplastic case is indistinguishable from the corresponding one in which the damage is ignored. As shown in Figure 15d, the evolution of damage in the metallic phase in the viscoplastic case is quite weak. Consequently, let the unidirectional metal/rubber-like composite be subjected to a uniaxial stress cyclic loading 0.5 ≤ F XX ≤ 1.5 at the off-axis angle φ = 10 • . The resulting response of the composite and the damage evolution in the fiber phase are shown in Figure 16. The graph shows that the damage increases rapidly and total failure of the fiber occurs after less than 3.5 cycles.

Application of the finite strain coupled viscoelastic-damage model
Let us consider a composite material that consists of a finite viscoelastic matrix whose free energy is given by Equation 44, which represents a single Maxwell element with the parameters given by Table 2. The damage parameters in Equation 53 are: D ∞ 0 = 1 and α dam = 1. The viscoelastic matrix is reinforced by continuous linearly elastic isotropic nylon fibers oriented in the 1-direction, whose Young's modulus and Poisson's ratio are 2GP a and 0.4, respectively. The volume fraction of the fibers is denoted by v f . Figure 17a shows the response of the monolithic finite viscoelastic material that is subjected to a uniaxial stress loading applied at a rate of 1s −1 for three values the damage parameter 1/η dam = 0, 1, 10. Figure 17b exhibits the corresponding damage evolutions in these cases. Referring to Equation 51, the value of 1/η dam = 0 corresponds to the case in which no damage takes place namely, D ≡ 0. The other parts of Figure 17 shows the response of the unidirectional composite and the corresponding damage evolution in the matrix phase for various amounts of fiber volume ratio and the above three values of damage parameter η dam . Changing the values of the latter parameters illustrates the effect of damage on the resulting behavior of the matrix and the overall response of the composite. The uniaxial stress loading is applied perpendicular to the fiber direction namely, in the transverse 2-direction. It should be noted that due to the large contrast between the moduli of the fiber and matrix, loading in the fiber direction would not exhibit the viscoelastic effects since in such a case the effect of the elastic fiber is dominant.  The responses and damage evolutions in Figure 17 were caused by uniaxial transverse stress loading applied at a rate of 1s −1 . Figure 18 exhibits the effect of applying the transverse loading at two different rates: 1s −1 and 0.01s −1 which is caused by the presence of the viscoelastic mechanism. In all cases shown in this figure, the damage parameter η dam = 0.1 is kept constant. Figure 18a,b show the response of the unreinforced finite viscoelastic matrix and the resulting damage evolution at these two rates, whereas the other portions of this figure exhibit the behavior of the nylon reinforced matrix with various     fiber volume ratios. It is interesting to observe that whereas the stress responses are sensitive to the rate of applied loading, the damage evolution in the matrix is not sensitive in the sense that up to the scale of the plot the behaviors caused by these two rates are indistinguishable.

Conclusions and Future Research
A finite strain micromechanical analysis has been presented which is capable of predicting the behavior of multiphase materials that are modeled, in the framework of continuum damage mechanics, by finite elastoplasticity, viscoplasticity and viscoelasticity coupled with damage. Applications were presented in various circumstances. The present applications can be readily extended to obtain the finite strain response of laminated composites that are subjected to in-plane loading.
The two-dimensional macroscopic Gurson's model which is suitable for the representation of porous materials, in which the porosity is oriented in the axial direction, was applied to both elastoplastic and viscoplastic phases and its predictions were compared with those provided by the HFGMC method. Extension to phases which are represented by the Gurson's three-dimensional model, in which the porosity is given by a spherical pore, is a subject for a future research. In the present research, Gurson's model was applied in its original form. The inclusion of damage threshold [12] according to which damage starts only at a critical value, and the incorporation of the improvements of Tvergaard [31]- [32] to obtain a closer agreement with numerical results of a periodic array of voids, and of Tvergaard and Needleman [33] to account for the effects of void nucleation and coalescence at failure is another subject for a future research.
In the present investigation, the finite strain constituent of the composite was modeled either as inelastic (time-dependent or time-independent) or viscoelastic constitutive relations. In Miehe and Keck [25] and Peric and Dettmer [30], finite strain generalized inelastic material models that combine elastic, inelastic and viscoelastic behavior was presented. Hence it is possible to generalize the present HFGMC model to include finite strain constituents in which the material behavior exhibits combined different rheological phenomena (elastic, elastoplastic, viscoplastic and viscoelastic). In the framework of infinitesimal strains, composites with phases that exhibit viscoelastic-viscoplastic behavior were investigated by Aboudi [34].
Due to the simplicity of the finite strain HFGMC model, it should not be difficult to link it to a finite element procedure in order to analyze composite structures (e.g., composite beams, plates and shells) undergoing large deformations. Indeed, the capability for such structural investigations has been already performed in the infinitesimal strain domain by Bednarcyk and Arnold [35] who presented a framework that enables coupled multiscale analysis of composite structures. To this end, they developed, free, finite element analysis-micromechanics analysis code (FEAMAC) software that couples the micromechanics analysis code of the generalized method of cells (MAC/GMC) with the commercial ABAQUS finite element software to perform micromechanics based finite element analysis such that the nonlinear composite material response at each integration point is modeled at each increment by MAC/GMC.
As to the small strain HFGMC method, It was recently coupled to ABAQUS software by Haj-Ali and Aboudi [36] to generate a nested local-global nonlinear finite element analysis of composite materials and structures. More recently, the hyperelastic HFGMC model has been coupled to the finite element ABAQUS software by Kim [37] to investigate the behavior of various types of tissue materials includ-ing the human arterial wall layers and porcine aortic valves leaflets. The results from this multiscale structural investigation were compared with the collagen fiber network (a model made of hyperelastic collagen and elastin layered finite elements) and Holzapfel et al. [38] and Gasser et al. [39] (hyperelastic anisotropic homogenized material) models. It was shown that the hyperelastic HFGMC is effective for the modeling of arteries especially when the collagen fiber network has a periodic microstructure.