Home | About | Journals | Submit | Contact Us | Français |

**|**Materials (Basel)**|**v.2(4); 2009 December**|**PMC5513566

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Finite Strain Coupled Elastoplasticity-Damage Models of Monolithic Materials
- 3. Finite Strain Coupled Viscoplasticity-Damage Models of Monolithic Materials
- 4. Finite Strain Coupled Viscoelasticity-Damage Model of Monolithic Materials
- 5. Finite Strain Micromechanical Analysis
- 6. Computational Procedures
- 7. Applications
- 8. Conclusions and Future Research
- References

Authors

Materials (Basel). 2009 December; 2(4): 1858–1894.

Published online 2009 November 18. doi: 10.3390/ma2041858

PMCID: PMC5513566

Faculty of Engineering, Tel Aviv University, Ramat Aviv 69978, Israel; E-Mail: li.ca.uat.gne@iduoba

Received 2009 October 14; Revised 2009 November 15; Accepted 2009 November 17.

Copyright © 2009 by the author.

Licensee Molecular Diversity Preservation International, Basel, Switzerland. This article is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license http://creativecommons.org/licenses/by/3.0/.

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.

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.

Let $\mathbf{X}$ and $\mathbf{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 $\mathbf{F}(\mathbf{X},t)$, $d\mathbf{x}=\mathbf{F}(\mathbf{X},t)d\mathbf{X}$. 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 $\mathbf{F}(\mathbf{X},t)$ in the form

$$\begin{array}{c}\hfill \mathbf{F}(\mathbf{X},t)={\mathbf{F}}^{e}(\mathbf{X},t){\mathbf{F}}^{p}(\mathbf{X},t)\end{array}$$

(1)

where ${\mathbf{F}}^{p}(\mathbf{X},t)$ and ${\mathbf{F}}^{e}(\mathbf{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

$$\begin{array}{c}\hfill \mathbf{C}={\mathbf{F}}^{T}\mathbf{F},\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\mathbf{C}}^{p}={\mathbf{F}}^{pT}{\mathbf{F}}^{p}\end{array}$$

(2)

where superscript *T* denotes the transpose operation. The left Cauchy-Green tensors $\mathbf{B}$ and ${\mathbf{B}}^{e}$ are defined by

$$\begin{array}{c}\hfill \mathbf{B}=\mathbf{F}{\mathbf{F}}^{T},\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\mathbf{B}}^{e}={\mathbf{F}}^{e}{\mathbf{F}}^{e}{}^{T}\end{array}$$

(3)

The logarithmic elastic strain ${\mathit{\u03f5}}^{e}$ is defined by

$$\begin{array}{c}\hfill {\mathit{\u03f5}}^{e}=\frac{1}{2}log\left[{\mathit{B}}^{e}\right]\end{array}$$

(4)

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.

The Hencky isotropic strain energy function per unit reference volume is defined by

$$\begin{array}{c}\hfill W=\frac{1}{2}{\mathit{\u03f5}}^{e}:\mathit{h}:{\mathit{\u03f5}}^{e}\end{array}$$

(5)

where the fourth-order tensor $\mathit{h}$ given by

$$\begin{array}{c}\hfill \mathit{h}=\lambda \mathit{I}\otimes \mathit{I}+2\mu {\mathit{I}}_{4}\end{array}$$

(6)

where *λ* and *μ* are material parameters, and $\mathit{I}$ and ${\mathit{I}}_{4}$ being the identity second and fourth-order tensors, respectively. The Kirchhoff stress * τ* is given by

$$\begin{array}{c}\hfill \mathit{\tau}=\frac{\partial W}{\partial {\mathit{\u03f5}}^{e}}=\mathit{h}:{\mathit{\u03f5}}^{e}\end{array}$$

(7)

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

$$\begin{array}{c}\hfill \Phi =\frac{1}{2}dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]:dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]-\frac{1}{3}\left[1+{D}^{2}-2Dcosh\left(\frac{\sqrt{3}({\tau}_{22}+{\tau}_{33})}{2{\sigma}_{y}}\right)\right]{\sigma}_{y}^{2}\end{array}$$

(8)

where $dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]$ is the deviator of Kirchhoff stress tensor, *D* denotes the porosity volume fraction (which represents the amount of damage) and ${\sigma}_{y}$ is the function that describes the hardening law of the material. For isotropic hardening it is given by

$$\begin{array}{c}\hfill {\sigma}_{y}={Y}_{0}+K\left(R\right)\end{array}$$

(9)

where ${Y}_{0}$ is the yield stress in simple tension and $K\left(R\right)$ describes the isotropic hardening law. For linear hardening: $K\left(R\right)=HR$.

The plastic flow rule is determined from Φ by employing the relation [12]:

$$\begin{array}{c}\hfill -\frac{1}{2}{L}_{v}\left[{\mathit{B}}^{e}\right]{\left[{\mathit{B}}^{e}\right]}^{-1}=\dot{\gamma}\frac{\partial \Phi}{\partial \mathit{\tau}}\end{array}$$

(10)

where ${L}_{v}$ is the Lie derivative of ${\mathit{B}}^{e}$ and *γ* is the consistency parameter. Here and in the following $\dot{Q}$ means the derivative of a quantity *Q* with respect to time or quasi time *t*. The Lie derivative of ${\mathit{B}}^{e}$ can be expressed in the form [13]

$$\begin{array}{c}\hfill {L}_{v}\left[{\mathit{B}}^{e}\right]=\mathit{F}{\dot{\mathit{C}}}^{p-1}{\mathit{F}}^{T}\end{array}$$

(11)

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 ${\left[{\mathit{B}}^{e}\right]}^{-1}$ on the left-hand-side replaced by ${\left(trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{B}}^{e}\right]\right)}^{-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]

$$\begin{array}{c}\hfill \dot{R}=-\frac{\dot{\gamma}}{1-D}\frac{\partial \Phi}{\partial {\sigma}_{Y}}\end{array}$$

(12)

whereas the evolution law for damage (porosity) is given by

$$\begin{array}{c}\hfill \dot{D}=\frac{2\dot{\gamma}}{\sqrt{3}}(D-{D}^{2}){\sigma}_{Y}sinh\left(\frac{\sqrt{3}({\tau}_{22}+{\tau}_{33})}{2{\sigma}_{y}}\right)\end{array}$$

(13)

Equations 10–13 together with the condition that $\Phi =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

$$\begin{array}{c}\hfill {\mathit{\u03f5}}_{n+1}^{e}={\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}-\Delta \gamma {\left.\frac{\partial \Phi}{\partial \mathit{\tau}}\right|}_{n+1}\end{array}$$

(14)

where ${\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ at the present $n+1$ increment is computed according to Equation 4 while employing ${\mathbf{B}}_{n}^{e}$ of the previous increment.

The fourth-order first tangent tensor $\mathit{R}$ that is needed in the micromechanical analysis is defined by

$$\begin{array}{c}\hfill \dot{\mathit{T}}=\mathit{R}\phantom{\rule{4pt}{0ex}}\dot{\mathit{F}}\end{array}$$

(15)

where $\mathit{T}$ is the first Piola-Kirchhoff stress tensor. It is related to the Kirchhoff stress tensor * τ* according to: $\mathit{T}={\mathit{F}}^{-1}\mathit{\tau}$. The first tangent tensor is determined according to

$$\begin{array}{ccc}\hfill \mathit{R}& =& \frac{\partial \mathit{T}}{\partial {\mathit{F}}_{n+1}}\hfill \\ & =& \frac{\partial {\mathit{F}}^{-1}}{\partial {\mathit{F}}_{n+1}}\mathit{\tau}+{\mathit{F}}^{-1}\frac{\partial \mathit{\tau}}{\partial {\mathit{F}}_{n+1}}\hfill \\ & =& \frac{\partial {\mathit{F}}^{-1}}{\partial {\mathit{F}}_{n+1}}\mathit{\tau}+{\mathit{F}}^{-1}\left[\frac{\partial \mathit{\tau}}{\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}:\frac{\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial {\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}:\frac{\partial {\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial {\mathit{F}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}:\frac{\partial {\mathit{F}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial {\mathit{F}}_{n+1}}\right]\hfill \end{array}$$

(16)

where

$$\begin{array}{ccc}\hfill {\left[\frac{\partial {\mathit{F}}^{-1}}{\partial \mathit{F}}\right]}_{ijkl}& =& -{F}_{ik}^{-1}{F}_{lj}^{-1}\hfill \\ \hfill \frac{\partial {\mathit{\u03f5}}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial {\mathit{B}}^{e\phantom{\rule{4pt}{0ex}}trial}}& =& \frac{1}{2}\frac{\partial log\left[{\mathit{B}}^{e\phantom{\rule{4pt}{0ex}}trial}\right]}{\partial {\mathit{B}}^{e\phantom{\rule{4pt}{0ex}}trial}}\hfill \\ \hfill {\left[\frac{\partial {\mathit{B}}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial {\mathit{F}}^{e\phantom{\rule{4pt}{0ex}}trial}}\right]}_{ijkl}& =& {\delta}_{ik}{F}_{jl}^{e\phantom{\rule{4pt}{0ex}}trial}+{\delta}_{jk}{F}_{il}^{e\phantom{\rule{4pt}{0ex}}trial}\hfill \\ \hfill {\left[\frac{\partial {\mathit{F}}^{e\phantom{\rule{4pt}{0ex}}trial}}{\partial \mathit{F}}\right]}_{ijkl}& =& {\delta}_{ik}{F}_{lj}^{p\phantom{\rule{4pt}{0ex}}-1}\hfill \end{array}$$

(17)

and ${\delta}_{ij}$ is the Kronecker delta. The tensor $\partial \mathit{\tau}/\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ 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 $\mathbf{R}$ (here and in all the following cases) is that the values of the first Piola-Kirchhoff $\mathit{T}$ that are obtained by the integration of Equation 15 must be identical to the values of $\mathit{T}={\mathit{F}}^{-1}\mathit{\tau}$ that are computed directly from Equation 7 in which $\mathbf{R}$ is not involved.

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):

$$\begin{array}{c}\hfill {\Psi}^{ed}=\frac{1}{2}{\mathit{\u03f5}}^{e}:(1-D)\mathit{h}:{\mathit{\u03f5}}^{e}\end{array}$$

(18)

where as before ${\mathit{\u03f5}}^{e}$ denotes the logarithmic elastic strain and *D* represents the damage variable such that $0\le D\le 1$. The Kirchhoff stress is determined according to, *cf.* Equation 7,

$$\begin{array}{c}\hfill \mathit{\tau}=\frac{\partial {\Psi}^{ed}}{\partial {\mathit{\u03f5}}^{e}}=(1-D)\mathit{h}:{\mathit{\u03f5}}^{e}\end{array}$$

(19)

The thermodynamical force *Y* which is conjugate to the damage variable *D* is given by

$$\begin{array}{c}\hfill Y=\frac{\partial {\Psi}^{ed}}{\partial D}=-\frac{1}{2}{\mathit{\u03f5}}^{e}:\mathit{h}:{\mathit{\u03f5}}^{e}\end{array}$$

(20)

By assuming isotropic hardening, the yield function Φ of this model is given by

$$\begin{array}{c}\hfill \Phi =\frac{1}{1-D}\sqrt{3{J}_{2}\left(\mathit{\tau}\right)}-{\sigma}_{y}\end{array}$$

(21)

where ${J}_{2}\left(\mathit{\tau}\right)=1/2\phantom{\rule{4pt}{0ex}}dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]:dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]=1/{2\phantom{\rule{4pt}{0ex}}\left|\right|dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]\left|\right|}^{2}$ being the second invariant of the Kirchhoff stress deviator tensor. This yield function readily determines the plastic flow rule which is given by Equation 10. The rate of hardening is given this time, *cf.* Equation 12, by

$$\begin{array}{c}\hfill \dot{R}=-\dot{\gamma}\frac{\partial \Phi}{\partial K}=\dot{\gamma}\end{array}$$

(22)

Finally, the dissipation function is taken in the form

$$\begin{array}{c}\hfill \Psi =\Phi +\frac{r}{(1-D)(s+1)}{\left(\frac{-Y}{r}\right)}^{s+1}\end{array}$$

(23)

where *r* and *s* are material constants. This function provides the evolution law of damage in the form

$$\begin{array}{c}\hfill \dot{D}=-\dot{\gamma}\frac{\partial \Psi}{\partial Y}=\frac{\dot{\gamma}}{1-D}{\left(\frac{-Y}{r}\right)}^{s}\end{array}$$

(24)

Equations 19–24 together with the condition that $\Phi =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 elastoplastic-damage model can be expressed in the form shown by Equation 15 where the fourth-order first tangent tensor $\mathbf{R}$ is determined by Equation 16. As in the Gurson’s model, the tensor $\partial \mathit{\tau}/\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ 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.

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 $\dot{\gamma}$ 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], $\dot{\gamma}$ in the coupled viscoplastic-damage Gurson’s model can be established in the form

$$\begin{array}{c}\hfill \dot{\gamma}=\frac{1}{\eta}\left[{\left(\frac{3\Phi}{a{\sigma}_{y}^{2}}+1\right)}^{1/2\u03f5}-1\right]H(\Phi )\end{array}$$

(25)

where the yield function Φ is given by Equation 8, $a=1+{D}^{2}-2Dcosh\left(\frac{\sqrt{3}({\tau}_{22}+{\tau}_{33})}{2{\sigma}_{y}}\right)$, *η* and *ϵ* are parameters that characterize the viscoplastic material, and $H(\Phi )$ is the Heaviside unit step function.

For Lemaitre coupled viscoplastic-damage model, $\dot{\gamma}$ is given by

$$\begin{array}{c}\hfill \dot{\gamma}=\frac{1}{\eta}\left[{\left(\frac{\Phi}{{\sigma}_{y}}+1\right)}^{1/\u03f5}-1\right]H(\Phi )\end{array}$$

(26)

where Φ is given by Equation 21. In the limit when $\eta \to 0$ and $\u03f5\to 0$, the viscoplastic model recovers the perfectly plastic von Mises model [12]. It should be noted that $\dot{\gamma}$ 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

$$\begin{array}{c}\hfill \dot{\gamma}=\frac{1}{\eta}{\left(\frac{\Phi}{{\sigma}_{y}}\right)}^{1/\u03f5}H(\Phi )\end{array}$$

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 $\mathbf{R}$ is given by Equation 16. In this latter equation, the instantaneous fourth-order tensor $\partial \mathit{\tau}/\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ is determined as follows. From Equation 14 we obtain that

$$\begin{array}{c}\hfill \Delta {\mathit{\u03f5}}^{e}={\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}-{\mathit{\u03f5}}_{n}^{e}-\Delta \gamma {\left.\frac{\partial \Phi}{\partial \mathit{\tau}}\right|}_{n+1}\end{array}$$

(27)

It follows from Equation 7 that for Gurson’s viscoplastic model

$$\begin{array}{c}\hfill \Delta \mathit{\tau}=\mathit{h}:\left(\Delta \mathit{\u03f5}-\Delta {\mathit{\u03f5}}^{vp}\right)\end{array}$$

(28)

where the following definitions are used

$$\begin{array}{c}\hfill \Delta \mathit{\u03f5}\equiv {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}-{\mathit{\u03f5}}_{n}^{e},\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\Delta {\mathit{\u03f5}}^{vp}\equiv \Delta \gamma {\left.\frac{\partial \Phi}{\partial \mathit{\tau}}\right|}_{n+1}\end{array}$$

(29)

Let us multiply $\mathit{h}:\Delta {\mathit{\u03f5}}^{vp}$ in Equation 28 by (*cf.* Boyd and Allen [28])

$$\begin{array}{c}\hfill \frac{\Delta {\mathit{\u03f5}}^{vp}:\mathit{h}:\Delta \mathit{\u03f5}}{\Delta {\mathit{\u03f5}}^{vp}:\mathit{h}:\Delta \mathit{\u03f5}}=1\end{array}$$

(30)

This readily yields that

$$\begin{array}{c}\hfill \Delta \mathit{\tau}=\left[\mathit{h}-\frac{\left(\mathbf{h}:\Delta {\mathit{\u03f5}}^{vp}\right)\otimes \left(\Delta {\mathit{\u03f5}}^{vp}:\mathbf{h}\right)}{\Delta {\mathit{\u03f5}}^{vp}:\mathit{h}:\Delta \mathit{\u03f5}}\right]:\Delta \mathit{\u03f5}\end{array}$$

(31)

It follows that the requested instantaneous tangent tensor for the viscoplastic Gurson’s model is of the form

$$\begin{array}{c}\hfill \frac{\partial \mathit{\tau}}{\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}=\mathit{h}-\frac{\left(\mathbf{h}:{\dot{\mathit{\u03f5}}}^{vp}\right)\otimes \left({\dot{\mathit{\u03f5}}}^{vp}:\mathbf{h}\right)}{{\dot{\mathit{\u03f5}}}^{vp}:\mathit{h}:\dot{\mathit{\u03f5}}}\end{array}$$

(32)

which can be employed in Equation 16.

For Lemaitre’s viscoplastic model, Equation 19 should be employed. Hence

$$\begin{array}{c}\hfill \Delta \mathit{\tau}=(1-D)\mathit{h}:\left(\Delta \mathit{\u03f5}-\Delta \gamma {\left.\frac{\partial \Phi}{\partial \mathit{\tau}}\right|}_{n+1}-\frac{\Delta D}{1-D}{\mathit{\u03f5}}_{n+1}^{e}\right)\end{array}$$

(33)

Let us denote

$$\begin{array}{c}\hfill \Delta {\mathit{\u03f5}}^{vpd}\equiv \Delta \gamma {\left.\frac{\partial \Phi}{\partial \mathit{\tau}}\right|}_{n+1}+\frac{\Delta D}{1-D}{\mathit{\u03f5}}_{n+1}^{e}\end{array}$$

(34)

Hence the instantaneous tangent tensor for this model is given by

$$\begin{array}{c}\hfill \frac{\partial \mathit{\tau}}{\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}}=(1-D)\mathit{h}-\frac{\left[(1-D)\mathbf{h}:{\dot{\mathit{\u03f5}}}^{vpd}\right]\otimes \left[{\dot{\mathit{\u03f5}}}^{vpd}:(1-D)\mathbf{h}\right]}{{\dot{\mathit{\u03f5}}}^{vpd}:(1-D)\mathit{h}:\dot{\mathit{\u03f5}}}\end{array}$$

(35)

which can be employed in Equation 16.

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

$$\begin{array}{c}\hfill \mathbf{F}(\mathbf{X},t)={\mathbf{F}}^{e}(\mathbf{X},t){\mathbf{F}}^{v}(\mathbf{X},t)\end{array}$$

(36)

where ${\mathbf{F}}^{v}$ is the counterpart of ${\mathbf{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:

$$\begin{array}{c}\hfill \psi ={\psi}^{EQ}+{\psi}^{NEQ}\equiv (1-D){\psi}_{0}^{EQ}+(1-D){\psi}_{0}^{NEQ}\end{array}$$

(37)

where ${\psi}_{0}^{EQ}$ and ${\psi}_{0}^{NEQ}$ are referred to as the effective free energy of the undamaged material, and *D* denotes the amount of damage such that $0\le D\le 1$.

The resulting Kirchhoff stresses are given by

$$\begin{array}{c}\hfill {\mathit{\tau}}^{EQ}=2\mathit{F}\frac{\partial {\psi}^{EQ}}{\partial \mathit{C}}{\mathit{F}}^{T}=(1-D){\mathit{\tau}}_{0}^{EQ}\end{array}$$

(38)

$$\begin{array}{c}\hfill {\mathit{\tau}}^{NEQ}=2{\mathit{F}}^{e}\frac{\partial {\psi}^{NEQ}}{\partial {\mathit{C}}^{e}}{\mathit{F}}^{e}{}^{T}=\frac{\partial {\psi}^{NEQ}}{\partial {\mathit{\u03f5}}^{e}}=(1-D){\mathit{\tau}}_{0}^{NEQ}\end{array}$$

(39)

where ${\mathit{C}}^{e}={\mathbf{F}}^{e}{}^{T}{\mathbf{F}}^{e}$ and ${\mathit{\tau}}_{0}^{EQ}$, ${\mathit{\tau}}_{0}^{NEQ}$ correspond to the effective Kirchhoff stresses of the undamaged material.

Let the left Cauchy-Green tensor $\mathit{B}$ be represented in terms of its eigenvalues:

$$\begin{array}{c}\hfill \mathit{B}=diag\phantom{\rule{4pt}{0ex}}[{b}_{1},{b}_{2},{b}_{3}]\end{array}$$

(40)

With $J=det\mathit{F}$, the volume preserving tensor $\overline{\mathit{B}}={J}^{-2/3}\mathit{B}$ can be accordingly represented in the form

$$\begin{array}{c}\hfill \overline{\mathit{B}}=diag\phantom{\rule{4pt}{0ex}}[{\overline{b}}_{1},{\overline{b}}_{2},{\overline{b}}_{3}]={\left({b}_{1}{b}_{2}{b}_{3}\right)}^{-1/3}\phantom{\rule{4pt}{0ex}}diag\phantom{\rule{4pt}{0ex}}[{b}_{1},{b}_{2},{b}_{3}]\end{array}$$

(41)

The finite strain elastic contribution can be modeled by the Ogden’s compressible material representation (Ogden [22], Holzapfel [21]) as follows

$$\begin{array}{c}\hfill {\psi}_{0}^{EQ}=\sum _{A=1}^{3}\omega \left({\overline{b}}_{A}\right)+\frac{{K}^{e}}{2}{\left(J-1\right)}^{2}\end{array}$$

(42)

where ${K}^{e}$ is the elastic bulk modulus and

$$\begin{array}{c}\hfill \omega \left({\overline{b}}_{A}\right)=\sum _{r=1}^{N}\frac{{\mu}_{r}^{e}}{{\alpha}_{r}^{e}}{\left({\overline{b}}_{A}\right)}^{{\alpha}_{r}^{e}/2}\end{array}$$

(43)

and ${\mu}_{r}^{e}$ and ${\alpha}_{r}^{e}$ are material parameters of the elastic element.

For Maxwell’s element, the free energy is represented by [19]:

$$\begin{array}{c}\hfill {\psi}_{0}^{NEQ}=\sum _{p=1}^{3}\frac{{\mu}_{p}^{v}}{{\alpha}_{p}^{v}}\left[{\left({\overline{b}}_{1}^{e}\right)}^{{\alpha}_{p}^{v}/2}+{\left({\overline{b}}_{2}^{e}\right)}^{{\alpha}_{p}^{v}/2}+{\left({\overline{b}}_{3}^{e}\right)}^{{\alpha}_{p}^{v}/2}-3\right]+\frac{{K}^{v}}{4}\left[{\left({J}^{e}\right)}^{2}-2log{J}^{e}-1\right]\end{array}$$

(44)

where

$$\begin{array}{c}\hfill {\mathit{B}}^{e}=diag\phantom{\rule{4pt}{0ex}}[{b}_{1}^{e},{b}_{2}^{e},{b}_{3}^{e}]\end{array}$$

(45)

and ${J}^{e}=\sqrt{{b}_{1}^{e}{b}_{2}^{e}{b}_{3}^{e}}$, ${\overline{b}}_{A}^{e}={\left({J}^{e}\right)}^{-2/3}{b}_{A}^{e}$, and ${\mu}_{p}^{v}$, ${\alpha}_{p}^{v}$, ${K}^{v}$ are material parameters. By employing Eqs. (37) and (39) the following expression for the principal values of ${\mathit{\tau}}_{0}^{NEQ}$ is obtained

$$\begin{array}{c}\hfill {\tau}_{0A}^{NEQ}=\sum _{p=1}^{3}{\mu}_{p}^{v}\left[\frac{2}{3}{\left({\overline{b}}_{A}^{e}\right)}^{{\alpha}_{p}^{v}/2}-\frac{1}{3}{\left({\overline{b}}_{B}^{e}\right)}^{{\alpha}_{p}^{v}/2}-\frac{1}{3}{\left({\overline{b}}_{C}^{e}\right)}^{{\alpha}_{p}^{v}/2}\right]+\frac{{K}^{v}}{2}\left[{\left({J}^{e}\right)}^{2}-1\right],\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}A,B,C=1,2,3\end{array}$$

(46)

The evolution equation in the present finite viscoelastic case that replaces Equation 10 is given by ([19])

$$\begin{array}{c}\hfill -\frac{1}{2}{L}_{v}\left[{\mathit{B}}^{e}\right]{\left[{\mathit{B}}^{e}\right]}^{-1}=\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\end{array}$$

(47)

where ${\eta}_{D}$ and ${\eta}_{V}$ are the deviatoric and volumetric viscosities, respectively, and presently (*cf.* Equation 11)

$$\begin{array}{c}\hfill {L}_{v}\left[{\mathit{B}}^{e}\right]=\mathit{F}{\dot{\mathit{C}}}^{v-1}{\mathit{F}}^{T}\end{array}$$

(48)

with ${\mathit{C}}^{v}={\mathit{F}}^{vT}{\mathit{F}}^{v}$. For elastic bulk behavior, $1/{\eta}_{V}=0$.

By employing the exponential mapping algorithm, Equation 47 is reduced to, *cf.* Equation 14,

$$\begin{array}{c}\hfill {\u03f5}_{n+1,A}^{e}={\u03f5}_{n+1,A}^{e\phantom{\rule{4pt}{0ex}}trial}-\Delta t{\left[\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{A}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\right]}_{n+1}\end{array}$$

(49)

where the principal values of the logarithmic strain ${\u03f5}_{A}^{e}$ are given by, *cf.* Equation 4, ${\u03f5}_{A}^{e}=1/2log\left({b}_{A}^{e}\right)$ and $\Delta t$ is the time increment between the current and previous step. Equation 49 forms a system of coupled nonlinear equations in the three unknowns: ${\u03f5}_{A}^{e}$, $A=1,2,3$. It readily provides, *cf.* Equation 27, that

$$\begin{array}{c}\hfill \Delta {\u03f5}_{A}^{e}={\u03f5}_{n+1,A}^{e\phantom{\rule{4pt}{0ex}}trial}-{\u03f5}_{n,A}^{e}-\Delta t{\left[\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{A}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\right]}_{n+1}\end{array}$$

(50)

The rate of damage evolution is given by Lin and Schomburg [24] and Miehe and Keck [25] in the form

$$\begin{array}{c}\hfill \dot{D}=\frac{\dot{z}}{{\eta}_{dam}}({D}^{\infty}-D)\end{array}$$

(51)

where

$$\begin{array}{c}\hfill \dot{z}=\sqrt{\frac{2}{3}}\left|\right|\dot{\mathit{H}}\left|\right|,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\mathit{H}=\frac{1}{2}log\mathit{C}\end{array}$$

(52)

with the saturation value

$$\begin{array}{c}\hfill {D}^{\infty}=\frac{1}{1+{D}_{0}^{\infty}exp\left(-{\beta}_{dam}/{\alpha}_{dam}\right)}\end{array}$$

(53)

and

$$\begin{array}{c}\hfill {\beta}_{dam}=\underset{0\le \xi \le t}{max}\sqrt{\frac{2}{3}}\left|\right|\mathit{H}\left(\xi \right)\left|\right|\end{array}$$

(54)

In these relations, ${\eta}_{dam}$, ${D}_{0}^{\infty}$ and ${\alpha}_{dam}$ are material parameters.

The fourth-order first tangent tensor ${\mathit{R}}^{NEQ}$ is determined by employing Equation 16. In this equation the expression $\partial {\mathit{\tau}}^{NEQ}/\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ is determined as follows. From Equation 39 the following expression can be established

$$\begin{array}{c}\hfill \Delta {\tau}_{A}^{NEQ}=(1-D)\Delta {\tau}_{0A}^{NEQ}-{\tau}_{0A}^{NEQ}\Delta D\end{array}$$

(55)

Let the second-order tensor $\mathit{M}$ be defined by

$$\begin{array}{c}\hfill \mathit{M}=\left[{M}_{AB}\right]\equiv \left[\frac{\partial {\tau}_{0A}^{NEQ}}{\partial {\u03f5}_{B}^{e}}\right],\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}A,B=1,2,3\end{array}$$

(56)

With

$$\begin{array}{c}\hfill {\tau}_{0A}^{NEQ}=dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{0A}^{NEQ}\right]+\frac{1}{3}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}_{0}^{NEQ}\right]\end{array}$$

(57)

the explicit components of $\mathit{M}$ are given by [19]:

$$\begin{array}{c}\hfill \frac{\partial dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{0A}^{NEQ}\right]}{\partial {\u03f5}_{A}^{e}}=\sum _{p=1}^{3}{\mu}_{p}^{v}{\alpha}_{p}^{v}\left[\frac{4}{9}{\left({\overline{b}}_{A}^{e}\right)}^{{\alpha}_{p}^{v}/2}+\frac{1}{9}{\left({\overline{b}}_{B}^{e}\right)}^{{\alpha}_{p}^{v}/2}+\frac{1}{9}{\left({\overline{b}}_{C}^{e}\right)}^{{\alpha}_{p}^{v}/2}\right]\end{array}$$

(58)

$$\begin{array}{c}\hfill \frac{\partial dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{0A}^{NEQ}\right]}{\partial {\u03f5}_{B}^{e}}=\sum _{p=1}^{3}{\mu}_{p}^{v}{\alpha}_{p}^{v}\left[-\frac{2}{9}{\left({\overline{b}}_{A}^{e}\right)}^{{\alpha}_{p}^{v}/2}-\frac{2}{9}{\left({\overline{b}}_{B}^{e}\right)}^{{\alpha}_{p}^{v}/2}+\frac{1}{9}{\left({\overline{b}}_{C}^{e}\right)}^{{\alpha}_{p}^{v}/2}\right]\end{array}$$

(59)

and

$$\begin{array}{c}\hfill \frac{\partial trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}_{0}^{NEQ}\right]}{\partial {\u03f5}_{A}^{e}}=3{K}^{v}{\left({J}^{e}\right)}^{2}\end{array}$$

(60)

In conjunction with Equation 50, we obtain that

$$\begin{array}{ccc}\hfill \Delta {\tau}_{A}^{NEQ}& =& (1-D){M}_{AB}\left\{{\u03f5}_{n+1,B}^{e\phantom{\rule{4pt}{0ex}}trial}-{\u03f5}_{n,B}^{e}-\Delta t\left[{\left.\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{B}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\right]}_{n+1}\right.\right\}\hfill \\ & -& {\tau}_{0A}^{NEQ}\Delta D\hfill \end{array}$$

(61)

Hence

$$\begin{array}{ccc}\hfill \Delta {\tau}_{A}^{NEQ}& =& (1-D){M}_{AB}\{{\u03f5}_{n+1,B}^{e\phantom{\rule{4pt}{0ex}}trial}-{\u03f5}_{n,B}^{e}-\Delta t{\left[\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{B}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\right]}_{n+1}\hfill \\ & -& {M}_{BQ}^{-1}{\tau}_{0Q}^{NEQ}\frac{\Delta D}{1-D}\}\hfill \end{array}$$

(62)

Let $\Delta {\u03f5}_{A}^{e}$ and $\Delta {\u03f5}_{A}^{ved}$ be defined by

$$\begin{array}{ccc}\hfill \Delta {\u03f5}_{A}& \equiv & {\u03f5}_{n+1,A}^{e\phantom{\rule{4pt}{0ex}}trial}-{\u03f5}_{n,A}^{e}\hfill \\ \hfill \Delta {\mathit{\u03f5}}_{A}^{ved}& \equiv & \Delta t{\left[\frac{1}{2{\eta}_{D}}dev\phantom{\rule{4pt}{0ex}}\left[{\tau}_{A}^{NEQ}\right]+\frac{1}{9{\eta}_{V}}trace\phantom{\rule{4pt}{0ex}}\left[{\mathit{\tau}}^{NEQ}\right]\right]}_{n+1}+{M}_{AB}^{-1}{\tau}_{0B}^{NEQ}\frac{\Delta D}{1-D}\hfill \end{array}$$

(63)

Therefore Equation 62 can be represented by

$$\begin{array}{c}\hfill \Delta {\tau}_{A}^{NEQ}=(1-D){M}_{AB}\left[\Delta {\u03f5}_{B}-\Delta {\u03f5}_{B}^{ved}\right]\end{array}$$

(64)

By following the previous derivation for the establishment of the tangent tensor in the viscoplastic case, the following expression can be obtained

$$\begin{array}{c}\hfill \frac{\partial {\tau}_{A}^{NEQ}}{\partial {\u03f5}_{B}^{e\phantom{\rule{4pt}{0ex}}trial}}=(1-D)\mathit{M}-\frac{\left[(1-D)\mathbf{M}:{\dot{\mathit{\u03f5}}}^{ved}\right]\otimes \left[{\dot{\mathit{\u03f5}}}^{ved}:(1-D)\mathbf{M}\right]}{{\dot{\mathit{\u03f5}}}^{ved}:(1-D)\mathit{M}:\dot{\mathit{\u03f5}}}\end{array}$$

(65)

With ${\tau}_{A}^{NEQ}$, $A=1,2,3$, and their derivatives given by Equations 39, 46 and 65, the derivative fourth-order tensor function $\partial {\mathit{\tau}}^{NEQ}/\partial {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ can be established by employing the explicit expressions given by [11]. Similar treatment applies to ${\mathit{\tau}}^{EQ}$ which leads to the establishment of $\partial {\mathit{\tau}}^{EQ}/\partial {\mathit{\u03f5}}_{n+1}^{e}$ and ${\mathit{R}}^{EQ}$. The first tangent tensor of the viscoelastic material is given by $\mathit{R}={\mathit{R}}^{NEQ}+{\mathit{R}}^{EQ}$. The special case of a hyperelastic material coupled with evolving damage can be obtained by disregarding the viscous effects so that ${\psi}^{NEQ}=0$.

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 $\mathbf{A}\left(\mathbf{Y}\right)$, where $\mathit{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 $\dot{\mathbf{F}}\left(\mathbf{Y}\right)$ at a material point $\mathbf{Y}$ within the repeating unit cell to the externally applied deformation gradient rate $\dot{\overline{\mathbf{F}}}$ in the form:

$$\dot{\mathbf{F}}\left(\mathbf{Y}\right)=\mathbf{A}\left(\mathbf{Y}\right):\dot{\overline{\mathbf{F}}}$$

(66)

It follows from Equation 66 and in conjunction with (15) that the local stress rate at this point is given by

$$\dot{\mathbf{T}}\left(\mathbf{Y}\right)=\mathbf{R}\left(\mathbf{Y}\right):\left(\mathbf{A}\left(\mathbf{Y}\right):\dot{\overline{\mathbf{F}}}\right)$$

(67)

Hence the resulting macroscopic constitutive equation for the multiphase composite undergoing large deformation is given by

$$\dot{\overline{\mathbf{T}}}={\mathbf{R}}^{*}:\dot{\overline{\mathbf{F}}}$$

(68)

where ${\mathbf{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 $\mathbf{R}\left(\mathbf{Y}\right)$ and the established deformation concentration tensor $\mathbf{A}\left(\mathbf{Y}\right)$ in the form

$${\mathbf{R}}^{*}=\frac{1}{{V}_{Y}}{\int}_{{V}_{Y}}\mathbf{R}\left(\mathbf{Y}\right)\phantom{\rule{4pt}{0ex}}\mathbf{A}\left(\mathbf{Y}\right)d{V}_{Y}$$

(69)

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 ${\mathbf{R}}^{*}$ of the composite is affected by the current value of damage variable through the instantaneous value of tangent tensors $\mathbf{R}\left(\mathbf{Y}\right)$ 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.

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

At time ${t}_{n}$, the local deformation gradient tensor ${\mathit{F}}_{n}$ and the elastic left Cauchy-Green tensor ${\mathit{B}}_{n}^{e}$ have been already established in the inelastic phase of the composite. In addition, ${\mathit{F}}_{n+1}$ is assumed to be known at time ${t}_{n+1}={t}_{n}+\Delta t$.

(1) Compute the local trial logarithmic strain ${\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ as follows:

$$\begin{array}{ccc}\hfill {\mathit{f}}_{n+1}& =& {\mathit{F}}_{n+1}{\mathit{F}}_{n}^{-1}\hfill \\ \hfill {\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}& =& {\mathit{f}}_{n+1}{\mathit{B}}_{n}^{e}{\mathit{f}}_{n+1}^{T}\hfill \\ \hfill {\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}& =& \frac{1}{2}log\left[{\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}\right]\hfill \end{array}$$

(70)

(2) With ${\mathit{\u03f5}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$, compute ${\mathit{\tau}}_{n+1}^{trial}$ 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, $\dot{\gamma}$ 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 $\mathit{R}$ of the constituent can be determined from Equation 16.

(4) With the established local value of $\mathit{R}$ in the inelastic phase, it is possible to proceed with the micromechanical analysis and compute the current concentration tensor $\mathit{A}$ that appears in Equation 66, see [1]. Consequently, Equation 69 can be employed to determine the instantaneous effective tangent tensor ${\mathit{R}}^{*}$ of the composite.

(5) The rate of the externally applied deformation gradient $\dot{\overline{\mathit{F}}}$ is determined in accordance with the prescribed type of loading. For a uniaxial deformation in the 1-direction, for example, ${\dot{\overline{F}}}_{11}$ is known while all other components of $\dot{\overline{\mathit{F}}}$ are zero. Hence, it is possible to compute the stress rate $\dot{\overline{\mathit{T}}}$ using Equation 68. In addition, by integrating the resulting global stress rate, $\overline{\mathit{T}}$ is obtained.

(6) Once $\dot{\overline{\mathit{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 ${\dot{\overline{F}}}_{11}$ only is known together with $\overline{\mathit{T}}$ the components of which are zero in the other directions, an iterative procedure is needed to determine the other components of $\dot{\overline{\mathit{F}}}$ from these conditions.

At time ${t}_{n}$, the local deformation gradient ${\mathit{F}}_{n}$ and the left Cauchy-Green deformation tensors ${\mathit{B}}_{n}^{e}$ have been already established, and ${\mathit{F}}_{n+1}$ is assumed to be known.

(1) From ${\mathit{F}}_{n}$ and ${\mathit{F}}_{n-1}$, the local right Cauchy-Green deformation tensors ${\mathit{C}}_{n}$ and ${\mathit{C}}_{n-1}$ can be readily determined. Hence $\dot{z}$ in Equation 52 can be determined, from which $\dot{D}$ that is given by Equation 51 can be integrated to provide the damage variable at time ${t}_{n+1}$:

$$\begin{array}{c}\hfill {D}_{n+1}=\frac{\Delta z{D}^{\infty}+{\eta}_{dam}{D}_{n}}{\Delta z+{\eta}_{dam}}\end{array}$$

(71)

Equation 54 yields

$${\beta}_{dam}^{n+1}=\left\{\begin{array}{cc}\hfill \sqrt{\frac{2}{3}}\left|\right|{\mathit{H}}_{n+1}\left|\right|,& if\phantom{\rule{4pt}{0ex}}\sqrt{\frac{2}{3}}\left|\right|{\mathit{H}}_{n+1}\left|\right|-{\beta}_{dam}^{n}>0\hfill \\ \hfill \sqrt{\frac{2}{3}}\left|\right|{\mathit{H}}_{n}\left|\right|,& \mathrm{otherwise}\hfill \end{array}\right.$$

(72)

(2) The local trial elastic left Cauchy-Green deformation tensor ${\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ can be computed from the second relation in Equation 70.

(3) From ${\mathit{B}}_{n+1}^{e\phantom{\rule{4pt}{0ex}}trial}$ the principal values ${b}_{A}^{e\phantom{\rule{4pt}{0ex}}trial}$ and the logarithmic strains ${\u03f5}_{A}^{e\phantom{\rule{4pt}{0ex}}trial}$ can be determined at time step $n+1$. The solution of the nonlinear equations (50) provides ${b}_{A}^{e}$, ${\u03f5}_{A}^{e}$ and ${\mathit{B}}^{e}$ at time step $n+1$.

(4) From the established ${b}_{A}^{e}$, the Kirchhoff stresses ${\tau}_{A}$ and tensor $\mathit{M}$ can be calculated from Equations 46 and 56, respectively.

(5) The local tangent tensors ${\mathit{R}}^{NEQ}$ and ${\mathit{R}}^{EQ}$ are determined from the procedure that leads to Equation 16 in conjunction with Equation 65.

(6) With the established local values of $\mathit{R}={\mathit{R}}^{NEQ}+{\mathit{R}}^{EQ}$ of the constituent, it is possible to proceed with the micromechanical analysis to compute the concentration tensor $\mathit{A}$ that appears in Equation 66, see [1]. Consequently, Equation 69 can be employed to determine the tangent tensor ${\mathit{R}}^{*}$.

(7) The rate of the externally applied deformation gradient $\dot{\overline{\mathit{F}}}$ is determined in accordance with the prescribed type of loading. For a uniaxial deformation in the 1-direction, for example, ${\dot{\overline{F}}}_{11}$ is known while all other components of $\dot{\overline{\mathit{F}}}$ are zero. Hence, it is possible to compute the stress rate $\dot{\overline{\mathit{T}}}$ using Equation 68. In addition, by integrating the resulting global stress rate, $\overline{\mathit{T}}$ is obtained.

(8) Once $\dot{\overline{\mathbf{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 ${\dot{\overline{F}}}_{11}$ only is known together with $\overline{\mathit{T}}$ the components of which are zero in the other directions, an iterative procedure is needed to determine the other components of $\dot{\overline{\mathit{F}}}$ from these conditions.

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 ${\eta}_{D}$ and $1/{\eta}_{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.
The parameters ${\mu}_{p}^{v}$ and ${\alpha}_{p}^{v}$, $p=1,2,2$ are the Ogden’s material constants, ${K}^{v}$ is the bulk modulus and ${\eta}_{D}$ is the viscoelastic constant of the viscous element with $1/{\eta}_{v}=0$.

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

$$\frac{3}{2}dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]:dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]=\left[1+{D}^{2}-2Dcosh\left(\frac{\sqrt{3}({\tau}_{22}+{\tau}_{33})}{2{Y}_{0}}\right)\right]{Y}_{0}^{2}$$

(73)

The von Mises yielding criterion is obtained from this equation by setting $D=0$. Suppose that in the Gurson’s model the initial yield surface ${\overline{T}}_{22}$ against ${\overline{T}}_{33}$ is requested, where $\overline{\mathit{T}}={\mathit{F}}^{-1}\mathit{\tau}$ being the first Piola-Kirchhoff stress tensor. Here all components ${\overline{T}}_{ij}$ are equal to zero except ${\overline{T}}_{22}$ and ${\overline{T}}_{33}$ that may take any combination. Hence let us represent these components in the form

$${\overline{T}}_{22}=Acos\theta ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\overline{T}}_{33}=Asin\theta $$

(74)

where *A* is the radial distance of a point located on the initial yield surface envelope ${\overline{T}}_{22}$-${\overline{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

$$\frac{A}{{Y}_{0}}=\sqrt{\frac{1+{D}^{2}-2Dcosh\left(\frac{\sqrt{3}A({F}_{22}cos\theta +{F}_{33}sin\theta )}{2{Y}_{0}}\right)}{{F}_{22}{cos}^{2}\theta +{F}_{33}{sin}^{2}\theta -{F}_{22}{F}_{33}sin\theta cos\theta}}$$

(75)

where in Equation 73, ${\tau}_{22}={F}_{22}{\overline{T}}_{22}$ and ${\tau}_{33}={F}_{33}{\overline{T}}_{33}$ 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}\approx 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 $\mathit{B}\left(\mathit{Y}\right)$ which relates the rate of the local Piola-Kirchhoff stress tensor $\dot{\mathit{T}}\left(\mathit{Y}\right)$ to the externally applied stress rate $\dot{\overline{\mathit{T}}}$, namely

$$\dot{\mathbf{T}}\left(\mathbf{Y}\right)=\mathbf{B}\left(\mathbf{Y}\right):\dot{\overline{\mathbf{T}}}$$

(76)

From Equations 66–68 this tensor can be shown to be given by

$$\mathbf{B}\left(\mathbf{Y}\right)=\mathbf{R}\left(\mathbf{Y}\right)\phantom{\rule{4pt}{0ex}}\mathbf{A}\left(\mathbf{Y}\right)\phantom{\rule{4pt}{0ex}}{\left[{\mathbf{R}}^{*}\right]}^{-1}$$

(77)

Corresponding to the above discussion, suppose that the initial yield envelope ${\overline{T}}_{22}$ against ${\overline{T}}_{33}$ is requested. Here too all the average (composite) Piola-Kirchhoff stress components $\overline{T}$ are equal to zero except ${\overline{T}}_{22}$ and ${\overline{T}}_{33}$ which can be represented by Equation 74. The initial yielding of any ductile phase of the HFGMC model is given by the von Mises criterion namely, it is given by Equation 73 in which $D=0$ is substituted. In addition, for initial yielding at very small strains, Equation 76 can be reduced to $\mathbf{T}=\mathbf{B}\overline{\mathbf{T}}$. Hence the initial yielding of a point whose polar angle is *θ* as predicted by the HFGMC is given by

$$\frac{A}{{Y}_{0}}=\underset{over\phantom{\rule{4pt}{0ex}}all\phantom{\rule{4pt}{0ex}}HFGMC\phantom{\rule{4pt}{0ex}}ductile\phantom{\rule{4pt}{0ex}}phases}{min}\sqrt{\frac{1}{{B}_{22}{cos}^{2}\theta +{B}_{33}{sin}^{2}\theta -{B}_{22}{B}_{33}sin\theta cos\theta}}$$

(78)

In Figure 1, comparisons between the initial yield surfaces ${\overline{T}}_{22}/{Y}_{0}$ against ${\overline{T}}_{33}/{Y}_{0}$ as predicted by the Gurson’s and HFGMC models are shown for three values of damage (porosity) values: $D=0.05,0.25$ and $0.4$. Also shown as a reference is the simple von Mises criterion of a homogeneous elastoplastic material namely: $\sqrt{\frac{3}{2}dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]:dev\phantom{\rule{4pt}{0ex}}\left[\mathit{\tau}\right]}={Y}_{0}$. It can be clearly observed that a fair correspondence between the two models exist. Next, Figure 2 shows the initial yield surfaces in the ${\overline{T}}_{11}$-${\overline{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 loading ${\overline{T}}_{11}$ and ${\overline{T}}_{22}={\overline{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 lines ${\overline{T}}_{11}-{\overline{T}}_{22}=\pm {Y}_{0}$ the first one of which passes the points: $({\overline{T}}_{11}/{Y}_{0}=1,{\overline{T}}_{22}/{Y}_{0}=0)$ and $({\overline{T}}_{11}/{Y}_{0}=0,{\overline{T}}_{22}/{Y}_{0}=-1)$, whereas the second one passes through the points: $({\overline{T}}_{11}/{Y}_{0}=-1,{\overline{T}}_{22}/{Y}_{0}=0)$ and $({\overline{T}}_{11}/{Y}_{0}=0,{\overline{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.

Comparisons between the initial yield surfaces ${\overline{T}}_{22}-{\overline{T}}_{33}$ as predicted by the Gurson’s (G) and HFGMC (HF) models. Also shown, as a reference, is the simple von Mises envelope in a homogeneous elastoplastic material (H). **...**

Comparisons between the initial yield surfaces ${\overline{T}}_{11}-{\overline{T}}_{22}$ as predicted by the Gurson’s (G) and HFGMC (HF) models. Also shown, as a reference, is the simple von Mises envelope in a homogeneous elastoplastic material (H). **...**

Comparisons between the initial yield surfaces ${\overline{T}}_{11}-{\overline{T}}_{22}={\overline{T}}_{33}$ as predicted by the Gurson’s (G) and HFGMC (HF) models. (a) $D=0.05$, (b) $D=0.25$, (c) $D=0.4$.

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 responses ${\overline{T}}_{11}-{\overline{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 responses ${\overline{T}}_{22}-{\overline{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 Figure 4 and Figure 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 Figure 4 and Figure 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}$. Figure 6 presents the uniaxial stress response ${\overline{T}}_{22}-{\overline{F}}_{22}$ to a transverse loading of the porous material as predicted by Gurson’s and HFGMC model. In the Gurson model, initial porosity of ${D}_{i}=0.05,0.25$ and $0.4$ evolve with loading to the final values: ${D}_{f}=0.084,0.35$ and $0.51$, respectively at which the applied deformation is terminated at ${\overline{F}}_{22}=1.5$. The corresponding HFGMC predictions are shown in this figure for these values of ${D}_{i}$ and ${D}_{f}$. It can be readily observed that Gurson’s predictions extend between HFGMC prediction for $D={D}_{i}$ and $D={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 response ${\overline{T}}_{22}-{\overline{F}}_{22}$ to a transverse loading of the composite whose porosity is *D* when: (1) ${D}_{m}=0$ (*i.e.*, no damage takes place in the matrix), (2) ${D}_{m}=D$ (*i.e.*, the damage in the matrix is kept equal to *D*) and, (3) ${D}_{m}$ evolves from ${D}_{mi}$ to the final value ${D}_{mf}$ when the loading reaches ${\overline{F}}_{22}=1.5$. In this last case and as indicated in Figure 7, for $D=0.05$, ${D}_{mi}=0.05$ reaching the final value of ${D}_{mf}=0.11$; for $D=0.25$, ${D}_{mi}=0.25$ reaching the final value of ${D}_{mf}=0.39$; and for $D=0.4$, ${D}_{mi}=0.4$ reaching the final value of ${D}_{mf}=0.54$. This figure well exhibits the effect of zero, constant and evolving damage in the ductile matrix of a porous material (composite) in which the amount of porosity is prescribed..

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.5MPa$, $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\equiv 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.

The rubber-like matrix is modeled as a hyperelastic compressible neo-Hookean material whose strain energy function is given by (Bonet and Wood [26])

$$W=\frac{\mu}{2}\left({I}_{1}-3-2logJ\right)+\frac{\lambda}{2}{(logJ)}^{2}$$

(79)

where ${I}_{1}$ is the first invariant of the right Cauchy-Green deformation tensor $\mathit{C}$, $J=det\mathit{F}$ and *λ* and *μ* are material constants. The second Piola-kirchhoff stress tensor $\mathit{S}$ is obtained from

$$\mathit{S}=2\frac{\partial W}{\partial C}$$

(80)

and the constitutive relations of this material can be established in the form given by Equation 15, where the instantaneous first tangent tensor $\mathit{R}$ of the material is given by

$$\mathit{R}=4\mathit{F}\frac{{\partial}^{2}W}{\partial \mathit{C}\partial \mathit{C}}{\mathit{F}}^{T}+\mathit{S}\otimes \mathit{I}$$

(81)

Figure 8(c) shows the uniaxial stress behavior of the monolithic hyperelastic material characterized by: $\lambda =980MPa$ and $\mu =30MPa$.

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 ${\mathbf{T}}^{\mathbf{X}}$, referred to the new coordinate system, are equal to zero except ${T}_{XX}$. In particular, $\varphi ={0}^{\circ}$ and ${90}^{\circ}$ 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\equiv 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 ($\varphi ={0}^{\circ}$), whereas loading in the transverse direction ($\varphi ={90}^{\circ}$) caused yielding at a later stage. Among the six off-axis angles at which the response of the composite is generated in the range of $1\le {F}_{XX}\le 3.5$, the highest yield stress is obtained at $\varphi ={60}^{\circ}$.

The response to uniaxial stress loading of a unidirectional metal/rubber-like composite in which no damage is taking place in the metallic phase ($D\equiv 0$) which is modeled by Lemaitre elastoplastic equations. The rotation *ϕ* around the 3-direction **...**

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\equiv 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 $\varphi ={0}^{\circ}$, the damage variable increases with applied loading until it reaches its final value of $D=1$. For the off-axis angle $\varphi ={20}^{\circ}$, 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 $\varphi ={0}^{\circ}$, ${10}^{\circ}$ and ${15}^{\circ}$, 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$.

The response to uniaxial stress loading at various off-axis angles *ϕ* of a unidirectional metal/rubber-like composite in the absence ($D=0$) and presence of evolving damage in the fiber phase which is modeled by Lemaitre elastoplastic equations. **...**

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 $\varphi ={20}^{\circ}$. It is interesting however to observe the response and damage evolution when the metallic/rubber-like composite is subjected to a cyclic loading such that $0.5\le {F}_{XX}\le 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 $\varphi ={25}^{\circ}$. 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.

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: $\eta =1s$ and $\u03f5=1$, respectively. The applied rate of loading is $1{s}^{-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 4. It should be noted that although the graphs of Figure 4 and Figure 12 exhibit similar behavior, the values of the responses are generally different. Indeed, for the axial loading in the 1-direction (along which the porosities are oriented) the predictions provided by the Gurson’s and HFGMC viscoplasticity are just like the time-independent elastoplastic case being quite close. However, for loading in the transverse 2-direction such closeness is not obtained. The responses in this case are shown in 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.

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: $\eta =1s$ and $\u03f5=1$, respectively. In all cases shown herein the loading was applied at a rate of $1{s}^{-1}$, and $r=4.5MPa$, $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\equiv 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.

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 uniaxial stress loadings ${F}_{XX}$ applied at off-axis angles $\varphi ={0}^{\circ}$ and ${10}^{\circ}$ 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 ($\varphi ={0}^{\circ}$) 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 $\varphi ={10}^{\circ}$ 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.

The response to uniaxial stress loading at two off-axis angles of a unidirectional metal/rubber-like composite in the absence ($D=0$) and presence of evolving damage in the fiber phase which is modeled by Lemaitre viscoplastic (VP) and elastoplastic (EP) **...**

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\le {F}_{XX}\le 1.5$ at the off-axis angle $\varphi ={10}^{\circ}$. 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.

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}^{\infty}=1$ and ${\alpha}_{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 $2GPa$ 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 $1{s}^{-1}$ for three values the damage parameter $1/{\eta}_{dam}=0,1,10$. Figure 17b exhibits the corresponding damage evolutions in these cases. Referring to Equation 51, the value of $1/{\eta}_{dam}=0$ corresponds to the case in which no damage takes place namely, $D\equiv 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 ${\eta}_{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 response to transverse uniaxial stress loading applied at a rate of $1{s}^{-1}$ of a composite that consists of a finite viscoelastic matrix reinforced by nylon elastic fibers oriented in the axial 1-direction. The figure shows the stress response **...**

The responses and damage evolutions in Figure 17 were caused by uniaxial transverse stress loading applied at a rate of $1{s}^{-1}$. Figure 18 exhibits the effect of applying the transverse loading at two different rates: $1{s}^{-1}$ and $0.01{s}^{-1}$ which is caused by the presence of the viscoelastic mechanism. In all cases shown in this figure, the damage parameter ${\eta}_{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.

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 including 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.

1. Aboudi J. Finite Strain Micromechanical Modeling of Multiphase Composites. Int. J. Multiscale Comput. Eng. 2008;6:411–434. doi: 10.1615/IntJMultCompEng.v6.i5.30. [Cross Ref]

2. Aboudi J. Finite Strain Micromechanical Analysis of Rubber-Like Matrix Composites Incorporating the Mullins Damage Effect. Int. J. Damage Mech. 2009;18:5–29. doi: 10.1177/1056789507081845. [Cross Ref]

3. Kachanov L.M. Introduction to Continuum Damage Mechanics. M. Nijhoff Publishres; Dordrecht, The Netherlands: 1986.

4. Lemaitre J., Chaboche J.L. Mechanics of Solid Materials. Cambridge University Press; Cambridge, UK: 1990.

5. Krajcinovic D. Damage Mechanics. Elsevier; New York, NY, USA: 1996.

6. Voyiadjis G.Z., Kattan P.I. Damage Mechanics. Taylor & Francis; Boca Raton, FL, USA: 2005.

7. Talreja R. A Continuum Mechanics Characterization of Damage in Composite Materials. Proc. R. Soc. Lond. 1985;A399:195–216. doi: 10.1098/rspa.1985.0055. [Cross Ref]

8. Gurson A.L. Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I—Yield Criteria and Flow Rules for Porous Ductile Media. J. Engng. Mater. Tech. 1977;99:2–15. doi: 10.1115/1.3443401. [Cross Ref]

9. Lemaitre J. A Continuous Damage Mechanics Model for Ductile Fracture. J. Engng. Mater. Tech. 1985;107:83–89. doi: 10.1115/1.3225775. [Cross Ref]

10. Lemaitre J. Coupled Elasto-Plasticity and Damage Constitutive Equations. Comp. Meth. Appl. Mech. Eng. 1985;51:31–49. doi: 10.1016/0045-7825(85)90026-X. [Cross Ref]

11. de Souza Neto E.A., Perić D., Owen D.R.J. Continuum Modelling and Numerical Simulation of Material Damage at Finite Strains. Arch. Comput. Meth. Engrg. 1998;5:311–384. doi: 10.1007/BF02905910. [Cross Ref]

12. de Souza Neto E.A., Perić D., Owen D.R.J. Computational Methods For Plasticity. Wiley; Chichester, UK:

13. Simo J.C., Hughes T.J.R. Computational Inelasticity. Springer-Verlag; New York, NY: 1998.

14. Weber G., Anand L. Finite Deformation Constitutive Equations and a Time Integration Procedure for Isotropic, Hyperelastic-Viscoplastic Solids. Comp. Meth. Appl. Mech. Eng. 1990;79:173–202. doi: 10.1016/0045-7825(90)90131-5. [Cross Ref]

15. Eterovic A.L., Bathe K.-J. A Hyperelastic Based Large Strain Elasto-Plastic Constitutive Formulation with Combined Isotropic-Kinematic Hardening Using the Logarithmic Stress and Strain Measures. Int. J. Num. Meth. Eng. 1990;30:1099–1114. doi: 10.1002/nme.1620300602. [Cross Ref]

16. Cuitino A., Ortiz M. A Material-Independent Method for Extending Stress Update Algorithms From Small-Strain Plasticity with Multiplicative Kinematics. Engng. Comp. 1992;9:437–451. doi: 10.1108/eb023876. [Cross Ref]

17. Simo J.C. Algorithms for Static and Dynamic Multiplicative Plasticity that Presuurve the Classical Return Mapping Schemes of the Infinitesimal Theory. Comp. Meth. Appl. Mech. Eng. 1992;99:61–112. doi: 10.1016/0045-7825(92)90123-2. [Cross Ref]

18. Perzyna P. Advances in Applied Mechanics. Vol. 9. Academic Press; New York, NY, USA: 1966. Fundamental Problems in Viscoplasticity; pp. 243–377.

19. Reese S., Govindjee S. A Theory of Finite Viscoelasticity and Numerical Aspects. Int. J. Solids Struct. 1998;35:3455–3482. doi: 10.1016/S0020-7683(97)00217-5. [Cross Ref]

20. Simo J.C. On a Fully Three-Dimensional Finite-Strain Viscoelastic Damage Model: Formulation and Computational Aspects. Comp. Meth. Appl. Mech. Eng. 1987;60:153–173. doi: 10.1016/0045-7825(87)90107-1. [Cross Ref]

21. Holzapfel G.A. Nonlinear Solid Mechanics. John Wiley; New York, NY, USA: 2000.

22. Ogden R.W. Non-linear Elastic Deformations. Ellis Horwood; Chichester, UK: 1984.

23. Govindjee S., Reese S. A Presentation and Comparison of Two Large Deformation Viscoelasticity Models. J. Engng. Mater. Tech. 1997;119:251–255. doi: 10.1115/1.2812252. [Cross Ref]

24. Lin R.C., Schomburg U. A Finite Elastic-Viscoelastic-Elastoplastic Material Law with Damage: Theoretical and Numerical Aspects. Comput. Methods Appl. Mech. Eng. 2003;192:1591–1627. doi: 10.1016/S0045-7825(02)00649-7. [Cross Ref]

25. Miehe C., Keck J. Superimposed Finite Elastic-Viscoelastic Stress Response with Damage in Filled Rubbery Polymers. Experiments, Modelling and Algorithmic Implementation. J. Mech. Phys. Solids. 2000;48:323–365. doi: 10.1016/S0022-5096(99)00017-4. [Cross Ref]

26. Bonet J., Wood R.D. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press; Cambridge, UK: 2008.

27. Aboudi J. Finite Strain Micromechanical Analysis for Thermoelastoplastic Multiphase Materials. J. Mech. Mater. Struct. 2008;3:809–829. doi: 10.2140/jomms.2008.3.809. [Cross Ref]

28. Boyd J.G., Allen D.H. High Temperature Constitutive Modeling-Theory and Applications. MD-vol. 26/AMD-vol. 121. ASME; New York, NY, USA: 1991. A Thermoviscoplastic Micromechanics Model for Short Fiber Composites; pp. 409–422.

29. Aboudi J., Pindera M.J. High-Fidelity Micromechanical Modeling of Continuously Reinforced Elastic Multiphase Materials Undergoing Finite Deformation. Math. Mech. of Solids. 2004;9:599–628. doi: 10.1177/1081286504038591. [Cross Ref]

30. Peric D., Dettmer W. A Computational Model for Generalized Inelastic Materials at Finite Strains Combining Elastic, Viscoelastic and Plastic Material Behaviour. Eng. Comp. 2003;20:768–787. doi: 10.1108/02644400310488862. [Cross Ref]

31. Tvergaard V. Infuence of Voids on Shear Band Instabilities under Plane Strain Conditions. Int. J. Fract. 1981;17:389–407. doi: 10.1007/BF00036191. [Cross Ref]

32. Tvergaard V. On Localization in Ductile Materials Containing Spherical Voids. Int. J. Fract. 1981;18:237–252.

33. Tvergaard V., Needleman A. Analysis of the Cup-Cone Fracture in a Round Tensile Bar. Acta Mettal. 1984;32:157–169. doi: 10.1016/0001-6160(84)90213-X. [Cross Ref]

34. Aboudi J. Micromechanically Established Constitutive Equations for Multiphase Materials with Viscoelasti-Viscoplastic Phases. Mech. Time-Dependent Mater. 2005;9:121–145. doi: 10.1007/s11043-005-1085-x. [Cross Ref]

35. Bednarcyk B.A., Arnold S.M. A Framework for Performing Multiscale Stochastic Progressive Failure Analysis of Composite Materials. NASA; Cleveland, Ohio, USA: 2007.

36. Haj-Ali R., Aboudi J. Nonlinear Micromechanical Formulation of the High Fidelity Generalized Method of Cells. Int. J. Solids Struct. 2009;46:2577–2592. doi: 10.1016/j.ijsolstr.2009.02.004. [Cross Ref]

37. Kim H.S. Ph.D. Dissertatation. Georgia Institute of Technology; Atlanta, GA, USA: 2009. Nonlinear Multiscale Anisotropic Material and Structural Models for Prosthetic and Native Aortic Heart Valves.

38. Holzapfel G.A., Gasser T.C., Ogden R.W. A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models. J. Elasticity. 2000;61:1–48. doi: 10.1023/A:1010835316564. [Cross Ref]

39. Gasser T.C., Ogden R.W., Holzapfel G.A. Hyperelastic Modelling of Arterial Layers with Distributed Collagen Fibre Orientations. J. Royal Soc. Interface. 2006;3:15–35. doi: 10.1098/rsif.2005.0073. [PMC free article] [PubMed] [Cross Ref]

Articles from Materials are provided here courtesy of **Multidisciplinary Digital Publishing Institute (MDPI)**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |