# FEM Implementation of the Coupled Elastoplastic/Damage Model: Failure Prediction of Fiber Reinforced Polymers (FRPs) Composites

I Ud Din, P. Hao, M. Aamir, Gérald Franz, Stéphane Panier

► To cite this version:

I Ud Din, P. Hao, M. Aamir, Gérald Franz, Stéphane Panier. FEM Implementation of the Coupled Elastoplastic/Damage Model: Failure Prediction of Fiber Reinforced Polymers (FRPs) Composites. Journal of Solid Mechanics, 2019, 11 (4), pp.842-853. 10.22034/jsm.2019.668617. hal-03845934

HAL Id: hal-03845934

https://u-picardie.hal.science/hal-03845934

Submitted on 9 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FEM Implementation of the Coupled Elastoplastic/Damage Model: Failure Prediction of Fiber Reinforced Polymers (FRPs) Composites

I. Ud Din<fnref n="1" />*, P. Hao<fnref n="1" />, M. Aamir<fnref n="2" />, G. Franz<fnref n="1" />, S. Panier<fnref n="1" />

Received 28 July 2019; accepted 1 October 2019

## ABSTRACT

The coupled damage/plasticity model for meso-level which is ply-level in case of Uni-Directional (UD) Fiber Reinforced Polymers (FRPs) is implemented. The mathematical formulations, particularly the plasticity part, are discussed in a comprehensive manner. The plastic potential is defined in effective stress space and the damage evolution is based on the theory of irreversible thermodynamics. The model which is illustrated here has been implemented by different authors previously but, the complete pre-requisite algorithm ingredients used in the implicit scheme implementation are not available in the literature. This leads to the complexity in its implementation. Furthermore, this model is not available as a built-in material constitutive law in the commercial Finite Element Method (FEM) softwares. To facilitate the implementation and understanding, all the mathematical formulations are presented in great detail. In addition, the elastoplastic consistent operator needed for implementation in the implicit solution scheme is also derived. The model is formularized in incremental form to be used in the Return Mapping Algorithm (RMA). The quasi-static load carrying capability and non-linearity caused by the collaborative effect of damage and plasticity is predicted with User MATerial (UMAT) subroutine which solves the FEM problem with implicit techniques in ABAQUS.

© 2019 IAU, Arak Branch. All rights reserved.

**Keywords:** Continuum damage mechanics; Plasticity/damage coupling; Meso-scale; Return mapping algorithm; Fiber reinforced polymers.

## 1 INTRODUCTION

COMPOSITES play an important role in the various industrial applications replacing the metals due to its better mechanical and structural functionalities. The mechanical behavior prediction of composites is very

<footnote n="1">¹Laboratoire des Technologies Innovantes, LTI-EA 3899, Université de Picardie Jules Verne, Amiens, France</footnote>

<footnote n="2">²School of Engineering, Edith Cowan University, Joondalup, Perth, Australia</footnote>

<footnote>*Corresponding author. Tel.: +33 758025599.</footnote>

<footnote>E-mail address: israr.uddin@u-picardie.fr (I. Ud Din).</footnote>
important in order to foresee the in-service performance of the structures. This virtual assessment of the mechanical behavior not only accelerates the design phase but also reduces the experimentation, and ultimately the cost of the product to the end user. However, anisotropic and heterogeneous nature of FRPs accompanied by complex failure and multi-damage modes has posed the mechanical behavior prediction a challenging task. Since the introduction of the advanced composites like FRPs, researchers have successfully developed mechanical behavior and failure prediction models based on the coupling of Continuum Damage Mechanics (CDM) and plasticity. These models possess some inherent advantages like the prediction of residual strength and stiffness as a function of load. The model developed by Ladevèze et al. [1] at the ply level is one of such models aimed to predict the mechanical response of continuous FRPs composites. Models that are based on the coupled damage/plasticity approach are capable to describe the non-linearity with higher precision that is due to the in-elastic strain and micro-damage induced in FRPs as compared to the elastic damage models. Since the nature of the composite structure is to accumulate damage continuously before the final collapse so the application of failure criteria alone does not seem sufficient until and unless accompanied by damage evolution as a continuous state variable [2,3]. Failure initiation and complete progression/evolution until final rupture shall be incorporated in the model in order to obtain an efficient, reliable and accurate prediction tool [4].

CDM based approach is well-suited for the intra-laminar failure and damage prediction due to the fact that it considers each micro-damage mode initiation/evolution and their interaction along with the corresponding degradation in the mechanical behavior. In other words, it conceptualizes the degradation in the elastic properties of a single lamina, intra-ply level, due to the formation and their interaction of the diffused micro-cracks and cohesion of fiber/matrix before the localization and conversion of the micro-cracks into a single macro-crack [5]. It also assumes the possible intermediate states for the damage corresponding to the magnitude of the loading and hence treating the damage as a continuous variable [6]. Employing damage as a continuous variable enables to predict the residual mechanical properties at a particular load level: for instance, residual strength or residual stiffness at a specific load level in quasi-static loading or corresponding to a load cycle in fatigue loading prior to the catastrophic collapse of the structure [7]. Ladevèze model [1] is of paramount importance in composite community for the failure analysis of FRPs due to its excellent prediction capability and robustness in terms of any arbitrary stacking sequence and orientation of the fibers for the quasi-static loading. In addition, this model has been employed successfully for impact and fatigue loading in UD FRPs as well as textile composites [5,8-11]. Nevertheless, this model is not available as a built-in material constitutive law in commercial FEM codes such as ABAQUS and ANSYS to the design engineers for the design and analysis of FRPs composites. The complete pre-requisite algorithm ingredient formulation required is not available in the literature that escalates the complexity in implementation specifically the plasticity part. Therefore, the present work is aimed to present the complete formulation along with the consistent operator required for the FEM implementation in the commercial code ABAQUS. After presenting the detailed formulation, the model is implemented in an implicit solver with user customized material behavior capability in ABAQUS known as UMAT.

## 2 LADEVEZE DAMAGE / PLASTICITY COUPLED MODEL

In CDM, all the micro-damages/voids are smeared-out or homogenized over the continuum and a new medium is assumed having the degraded stiffness. This methodology was first proposed by Kachanov investigating the creep analysis of isotropic materials [12]. It is based on the theory of irreversible thermodynamic process considering internal state variables. These internal state variables can be expressed by a set of scalar or in tensorial notion with their conjugate driving force variables [13]. The choice of the state variables depends on the material damage modes. These memorize the history that what happened to the material throughout the loading range [7].

### 2.1 Plane stress anisotropic damage constitutive law

FRPs based composites fail in a variety of mechanisms which ranging from micro-level to macro-level. The common micro-level failure mechanisms comprise of fiber fracture, fiber micro-buckling (kinking) in compression, fiber/matrix debonding, matrix micro-cracking [14]. These micro-level failures in the layers then act as a foundation for the subsequent meso-level (e.g. transverse matrix cracking in a ply) and macro-level failure on laminate level such as delamination [5]. The strain energy of the damage orthotropic material denoted by $E_D$ can be expressed based on the plane stress assumption as [8,1]:
$$E_D = \frac{1}{2} \left[ \frac{\sigma_{11}^2}{E_{11}^0 (1-d_1)} - \frac{2\nu_{12}\sigma_{11}\sigma_{22}}{E_{11}^0 (1-d_1)} + \frac{\langle\sigma_{22}\rangle_+^2}{E_{22}^0 (1-d_2)} + \frac{\langle\sigma_{22}\rangle_-^2}{E_{22}^0} + \frac{(\tau_{12})^2}{G_{12}^0 (1-d_{12})} \right] \quad (1)$$

where $E_{11}^0, E_{22}^0, G_{12}^0, \nu_{12}$ denote the virgin material (undamaged) in-plane elastic constants. No degradation is assumed in the Poisson ratios $\nu_{12}$. And, $d_1, d_2, d_{12}$ are the scalar damage state variables in fiber direction, transverse direction and in-plane shear that degrade the virgin material elastic constants in the respective directions proportional to the state of stress. Based on the plane stress assumption the components of stress tensor are denoted as $\sigma_{11}, \sigma_{22}$ and $\tau_{12}$. In case of transverse compression, it is assumed that cracks are closed. Therefore, crack closure can be defined as:

$$\langle a \rangle_+ = a \quad \text{if } a \ge 0; \text{ otherwise } \langle a \rangle_+ = 0 \qquad (2)$$

$$\langle a \rangle_- = a \quad \text{if } a \le 0; \text{ otherwise } \langle a \rangle_- = 0 \qquad (2)$$

The damage orthotropic elastic constitutive law is thus obtained:

$$\boldsymbol{\epsilon}^e = \frac{\partial E_D}{\partial \boldsymbol{\sigma}} = \mathbf{S}^e : \boldsymbol{\sigma} \Leftrightarrow \left\{ \begin{array}{l} \epsilon_{11}^e = \frac{\sigma_{11}}{E_{11}^0 (1-d_1)} - \frac{\nu_{12}}{E_{11}^0 (1-d_1)} \sigma_{22} \\ \epsilon_{22}^e = \frac{\langle\sigma_{22}\rangle_+}{E_{22}^0 (1-d_2)} + \frac{\langle\sigma_{22}\rangle_-}{E_{22}^0} - \frac{\nu_{12}}{E_{11}^0} \sigma_{11} \\ \gamma_{12}^e = \frac{\tau_{12}}{G_{12}^0 (1-d_{12})} \end{array} \right\} \qquad (3)$$

where $\boldsymbol{\sigma}$ and $\boldsymbol{\epsilon}^e$ are stress and elastic strain tensors respectively which are treated as vectors in Voigt notation due to symmetry. $\mathbf{S}^e$ is the compliance matrix for the case of orthotropic material in plane stress case is presented in terms of elastic constants as:

$$\mathbf{S}^e = (\mathbf{D}^e)^{-1} = \begin{bmatrix} C_{11} & C_{12} & 0 \\ C_{12} & C_{22} & 0 \\ 0 & 0 & C_{33} \end{bmatrix} = \begin{bmatrix} \frac{1}{E_{11}^0} & -\frac{\nu_{12}}{E_{11}^0} & 0 \\ -\frac{\nu_{21}}{E_{22}^0} & \frac{1}{E_{22}^0} & 0 \\ 0 & 0 & \frac{1}{G_{12}^0} \end{bmatrix} \qquad (4)$$

The associated conjugate thermodynamic force variables denoted by $Y_2$ and $Y_{12}$ of the anisotropic damage $d_2$ and $d_{12}$ are determined from the Helmholtz free energy density $\Psi$. These thermodynamic forces are basically the driving forces which accumulate the damage as the external loading on the structure is increased.

$$Y_2 = \rho \left. \frac{\partial \Psi}{\partial d_2} \right|_{\sigma, d_{12}} = \rho \left. \frac{\partial E_D}{\partial d_2} \right|_{\sigma, d_{12}} = \frac{\sigma_{22}^2}{2E_{22}^0 (1-d_2)^2} \qquad (5)$$

$$Y_{12} = \rho \left. \frac{\partial \Psi}{\partial d_{12}} \right|_{\sigma, d_2} = \rho \left. \frac{\partial E_D}{\partial d_{12}} \right|_{\sigma, d_2} = \frac{\tau_{12}^2}{2G_{12}^0 (1-d_{12})^2} \qquad (6)$$

A linear coupling of the two conjugate forces is necessary to obtain an equivalent thermodynamic force $\mathbf{Y}(t)$ to be used in the damage accumulation [1]:
$$ \underline{Y}(t) = \max_{s \le t} (\sqrt{Y_{12}(s)} + bY_{2}(s)) \qquad (7) $$

$$ \underline{Y}_{-2}(t) = \max_{s \le t} (\sqrt{\underline{Y}_{2}(s)}) \qquad (8) $$

Here b is the shear and transverse damage coupling parameter and is determined experimentally. The damage evolution laws for the shear and transverse loading are determined by fitting the data points from the hysteresis loading-unloading stress-strain curve knowing the damage variables at each cycle and the corresponding stress level in pure shear laminate [($\pm 45^\circ$)]$_{2s}$ and angle-ply [($\pm 67.5^\circ$)]$_{2s}$ as proposed in [1]. A linear damage evolution law can be obtained as given below:

$$ d_{2} = \frac{\langle \underline{Y}(t) - Y_{2(0)} \rangle_{+}}{Y_{2(c)}} \quad \text{if } d_{2} < 1, \underline{Y}_{2}(t) \le Y_{s} \qquad (9) $$

$$ d_{12} = \frac{\langle \underline{Y}(t) - Y_{12(0)} \rangle_{+}}{Y_{12(c)}} \quad \text{if } d_{12} < 1, \underline{Y}_{2}(t) \le Y_{s} \qquad (10) $$

It shall be noted that the damage evolution in $d_1$ is accomplished with the help of Eq.(25) proposed in [8]. Though for the sake of completeness, the model parameters are briefly discussed here but the detailed description, physical meaning and experimental identification in the current formulations can be found in [1,8,15]. $Y_{2(0)}$ and $Y_{12(0)}$ are the threshold values of the thermodynamic forces at which intra-laminar micro-damage commences in the transverse and in-plane shear loading respectively. Moreover, $Y_{12(c)}$, $Y_s$ are the critical values of the thermodynamic forces in the shear and transverse direction, respectively, where the damage accumulation achieves its saturation. So, the material is regarded as damaged completely and it cannot sustain additional loading beyond the damage saturation values.

The effective stress tensor $\tilde{\sigma}$ can be found from the nominal Cauchy stress tensor $\sigma$ using the damage effect tensor $M$ that in case of reinforced materials with stiff and strong fibers can be accurately represented by a second order tensor in which the principal directions are aligned with the material direction:

$$ \tilde{\sigma} = M : \sigma \qquad (11) $$

$$ \text{Here } M = \text{diag} \left[ \frac{1}{1-d_1}, \frac{1}{1-d_2}, \frac{1}{1-d_{12}} \right] $$

## 2.2 Plasticity modeling

In case of the coupled plasticity/damage model, the plastic potential shall be expressed in the effective stress space [16]. Therefore, the plasticity potential used by Ladevèze et al., Eq.(12), is rendered in the quadratic form as given in Eq.(13). It is noteworthy that the quadratic form is easier and compact to implement with respect to the algorithm ingredients formulation derivation as extensively discussed in [17] for the generalized elastoplastic models with some particularizations for special cases. The elastoplastic formulation in [17] was extended to elastoplastic coupled damage model by the authors in [18] for a novel coupled elastoplastic damage model for FRPs composites. In the present work, a similar approach for the coupled damage/plasticity is followed as illustrated in [18]. The quadratic formulation as presented here for the damage/plasticity coupled model is very efficient to tailor to other plasticity functions and damage effect tensors accordingly by only changing the mapping matrices involved. This research work is intended with aim to facilitate the implementation of Ladevèze model developed in [1].

$$ f(\tilde{\sigma}, \tilde{p}) = \sqrt{(\tilde{\sigma}_{12}^2 + a^2 \tilde{\sigma}_{22}^2) - \bar{\sigma}_y}(\tilde{p}) = 0 \qquad (12) $$
$$F(\tilde{\sigma}, \tilde{p}) = \frac{1}{2} \tilde{\sigma}^T : \mathbb{P} : \tilde{\sigma} - \bar{\sigma}_y^2(\tilde{p}) \qquad (13)$$

Here $\bar{\sigma}_y(\tilde{p})$ is the evolving yield limit which increases as the plastic strain $\tilde{p}$ increases with load in case of strain hardening, $\bar{\sigma}_0$ is the initial yield limit, $\beta$ is the hardening law coefficient and $\alpha$ is the isotropic hardening law exponent. They are correlated here as the power law i.e $\bar{\sigma}_y(\tilde{p}) = \beta \tilde{p}^{\alpha} + \bar{\sigma}_0$. The parameter that establishes coupling between the transverse plasticity and in-plane shear plasticity is denoted by $a$. The mapping matrix that depends only on the constant ($a$) of the plasticity potential of the model is noted by $\mathbb{P}$ and can be read as:

$$\mathbb{P} = 2 \begin{bmatrix} 0 & 0 & 0 \\ 0 & a^2 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

### 2.2.1 Return mapping algorithm (RMA)

For the implementation of the model in FEM, the RMA [19] is followed which is a standard iterative problem to find the shortest distance between a point (trial state) and a convex set specified by the plasticity potential Eq. (13). The necessary steps for the implementation in ABAQUS UMAT are elaborated here.

**Step 1:**

Initial conditions at the start of the FEM increment:

$$\Delta\boldsymbol{\epsilon}, \boldsymbol{\epsilon}_{(n)}^e, \tilde{\boldsymbol{\epsilon}}_{(n)}^p, \tilde{p}_{(n)}, \boldsymbol{\sigma}_{(n)}, d_{1,n}, d_{2,n}, d_{12,n}$$

**Step 2:**

Trial stress based on the elastic constitutive law:

$$\boldsymbol{\epsilon}_{(n+1)}^{\text{tr}} = \boldsymbol{\epsilon}_{(n)}^{e} + \Delta\boldsymbol{\epsilon} \\ \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\text{tr}} = \mathbb{D}^{e}: \boldsymbol{\epsilon}_{(n+1)}^{\text{tr}} = \tilde{\boldsymbol{\sigma}}_{(n)} + \mathbb{D}^{e}: \Delta\boldsymbol{\epsilon}$$

$\mathbb{D}^e$ is the elastic stiffness tensor for the transversely isotropic UD ply:

$$\mathbb{D}^e = \begin{bmatrix} D_{11} & D_{12} & 0 \\ D_{12} & D_{22} & 0 \\ 0 & 0 & D_{33} \end{bmatrix} = \begin{bmatrix} \frac{E_{11}^0}{(1-\nu_{12}\nu_{21})} & \frac{\nu_{12}E_{22}^0}{(1-\nu_{12}\nu_{21})} & 0 \\ \frac{\nu_{12}E_{22}^0}{(1-\nu_{12}\nu_{21})} & \frac{E_{22}^0}{(1-\nu_{12}\nu_{21})} & 0 \\ 0 & 0 & G_{12}^0 \end{bmatrix} \qquad (14)$$

**Step 3:**

If $F(\tilde{\boldsymbol{\sigma}}_{(n+1)}, \tilde{p}_{(n)}) \le 0$ then the FEM increment is elastic and update the following quantities in UMAT:

$$\tilde{\boldsymbol{\sigma}}_{(n+1)} = \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\text{tr}}$$

$$\boldsymbol{\epsilon}_{(n+1)} = \boldsymbol{\epsilon}_{(n+1)}^{\text{tr}}$$

$$\tilde{\boldsymbol{\epsilon}}_{(n+1)}^{\text{p}} = \tilde{\boldsymbol{\epsilon}}_{(n)}^{\text{p}}$$

$$\tilde{p}_{(n+1)} = \tilde{p}_{(n)}$$
Otherwise, the current loading step is plastic and the plasticity potential is now based on the corrected stress tensor which is unknown for the moment:

$$F(\tilde{\sigma}_{(n+1)}, \tilde{p}_{(n+1)}) = \frac{1}{2}(\tilde{\sigma}_{(n+1)})^T : \mathbb{P} : \tilde{\sigma}_{(n+1)} - \bar{\sigma}_y^2 (\tilde{p}_{(n+1)}) \quad (15)$$

Step 4:

Find the Lagrange multiplicative using Newton-Raphson iterative scheme. It is required to derive the returned stress tensor (corrected stress tensor) in terms of the trial stress tensor. Therefore based on the elastic constitutive law we get:

$$\begin{aligned} \tilde{\sigma}_{(n+1)} &= \mathbb{D}^e : \boldsymbol{\varepsilon}_{(n+1)}^e = \mathbb{D}^e : (\boldsymbol{\varepsilon}_{(n+1)}^{\text{tr}} - \mathbf{d}\boldsymbol{\varepsilon}_{(n+1)}^p) \\ \tilde{\sigma}_{(n+1)} &= \tilde{\sigma}_{(n+1)}^{\text{tr}} - \mathbb{D}^e : \mathbf{d}\boldsymbol{\varepsilon}_{(n+1)}^p \end{aligned} \quad (16)$$

$\mathbf{d}\boldsymbol{\varepsilon}_{(n+1)}^p$ is the differential in plastic strain vector and is given by: $\mathbf{d}\boldsymbol{\varepsilon}_{(n+1)}^p = \mathbb{M} : \mathbf{d}\tilde{\boldsymbol{\varepsilon}}_{(n+1)}^p$ where $\mathbb{M}$ is the damage effect tensor defined previously in Eq.(11). The effective plastic strain increment is determined by the flow rule i.e $\mathbf{d}\tilde{\boldsymbol{\varepsilon}}^p = \Delta\lambda \frac{\partial F}{\partial \tilde{\boldsymbol{\varepsilon}}}$ where $\Delta\lambda$ is the Lagrange's plastic multiplicative. It is important to note that when the plasticity function in Eq.(12) is used, then $\Delta\lambda = d\tilde{\boldsymbol{p}}$. However, in case of quadratic plasticity function, these two quantities are not equal. Therefore, equating the plastic work to the equivalent plastic work, the expression for accumulated plastic strain can be obtained. Eq.(17) describes the accumulated plastic strain and its matrix notation counterpart.

$$d\tilde{\boldsymbol{p}} = \sqrt{\left(d\tilde{\gamma}_{12}\right)^2 + \frac{1}{a^2}\left(d\tilde{\varepsilon}_{22}^p\right)^2} = \sqrt{\left(\mathbf{d}\tilde{\boldsymbol{\varepsilon}}^p\right)^T : \mathbb{Z} : \mathbf{d}\tilde{\boldsymbol{\varepsilon}}^p} \quad (17)$$

It shall be noted that $d\tilde{\gamma}_{12}$ is the differential in the engineering shear strain ($d\tilde{\gamma}_{12} = 2d\tilde{\varepsilon}_{12}^p$) and $\mathbb{Z}$ is the mapping matrix referred below:

$$\mathbb{Z} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & \frac{1}{a^2} & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

The normal vector $\mathbf{n}_F$ to the scalar-valued plastic domain having the same principal axes to that of $\tilde{\boldsymbol{\sigma}}$ is defined as: $\mathbf{n}_F = \frac{\partial F}{\partial \tilde{\boldsymbol{\sigma}}} = \mathbb{P} \cdot \tilde{\boldsymbol{\sigma}}$. The effective plastic strain is then read as: $\mathbf{d}\tilde{\boldsymbol{\varepsilon}}^p = \Delta\lambda \mathbb{P} \cdot \tilde{\boldsymbol{\sigma}}$. Eq.(17) can be represented in terms of the returned stress tensor:

$$d\tilde{\boldsymbol{p}} = \Delta\lambda \sqrt{(\mathbb{P} \cdot \tilde{\boldsymbol{\sigma}})^T : \mathbb{Z} : (\mathbb{P} \cdot \tilde{\boldsymbol{\sigma}})} \quad (18)$$

Put the above expression in Eq.(16), the returned stress tensor can be computed in terms of the trial stress tensor:

$$\begin{aligned} \tilde{\boldsymbol{\sigma}}_{(n+1)} &= \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\text{tr}} - \Delta\lambda \mathbb{D}^e \cdot \mathbb{M} \cdot \mathbb{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)} \\ \tilde{\boldsymbol{\sigma}}_{(n+1)} &= (\mathbb{I} + \Delta\lambda \mathbb{D}^e \cdot \mathbb{M} \cdot \mathbb{P})^{-1} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\text{tr}} \end{aligned} \quad (19)$$

Here $\mathbb{I}$ is the identity matrix. Eq.(15) finally renders as:
$$F(\Delta\lambda) = \frac{1}{2} \left[ (\mathbf{I} + \Delta\lambda \mathbf{D}^e \cdot \mathbf{M} \cdot \mathbf{P})^{-1} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\mathrm{tr}} \right]^T \mathbf{P} \left[ (\mathbf{I} + \Delta\lambda \mathbf{D}^e \cdot \mathbf{M} \cdot \mathbf{P})^{-1} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)}^{\mathrm{tr}} \right] - \bar{\sigma}_y^2 (\tilde{p}_{(n+1)}) \quad (20)$$

The accumulated plastic strain in the above expression is $\tilde{p}_{(n+1)} = \tilde{p}_{(n)} + \Delta\lambda \sqrt{(\mathbf{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)})^T \mathbf{Z} \cdot (\mathbf{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)})}$. As the only unknown in Eq.(20) is the Lagrange's plastic multiplicator $\Delta\lambda$ which can be solved using Newton-Raphson iterative scheme for finding the approximate roots of the non-linear equation.

### 2.2.1.1 Newton-raphson iterative scheme

The $k^{th}$ guess for the solution of $\Delta\lambda$ is obtained as: $\Delta\lambda^{(k)} = \Delta\lambda^{(k-1)} + d \Delta\lambda^{(k)}$ where $d \Delta\lambda^{(k)}$ is given as follows:

$$d \Delta\lambda^{(k)} = - \left. \frac{F(\Delta\lambda)}{dF/d\Delta\lambda} \right|^{k-1} \quad (21)$$

The differentiation of Eq.(20) with respect to $\Delta\lambda$ is needed for the Newton-Raphson iterative scheme in Eq.(21). Following are the complete derivatives needed for UMAT:

$$\frac{dF}{d\Delta\lambda} = (\tilde{\boldsymbol{\sigma}}_{(n+1)})^T \mathbf{P} \frac{d\tilde{\boldsymbol{\sigma}}_{(n+1)}}{d\Delta\lambda} - 2\bar{\sigma}_y (\tilde{p}_{(n+1)}) \frac{d\bar{\sigma}_y}{d\Delta\lambda} (\tilde{p}_{(n+1)}) \quad (22)$$

$$\text{where } \frac{d\bar{\sigma}_y}{d\Delta\lambda} (\tilde{p}_{(n+1)}) = \alpha\beta \tilde{p}_{(n+1)}^{\alpha-1} \frac{d\tilde{p}_{(n+1)}}{d\Delta\lambda}$$

$$\frac{d\tilde{p}_{(n+1)}}{d\Delta\lambda} = \frac{(\mathbf{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)})^T \mathbf{Z} \mathbf{P} \left( \tilde{\boldsymbol{\sigma}}_{(n+1)} + \Delta\lambda \frac{d\tilde{\boldsymbol{\sigma}}_{(n+1)}}{d\Delta\lambda} \right)}{\sqrt{(\mathbf{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)})^T \mathbf{Z} \cdot (\mathbf{P} \cdot \tilde{\boldsymbol{\sigma}}_{(n+1)})}} \quad (23)$$

### 2.2.1.2 Algorithm consistent operator

For implementation of the above model in ABAQUS UMAT, the elastoplastic consistent tangent operator is required to determine in order to achieve convergence. Therefore, following [17], the elastoplastic consistent tangent matrix denoted by $\mathbf{D}^{ep}$ (Consistent Jacobian: $\frac{\partial\Delta\tilde{\boldsymbol{\sigma}}_{(n+1)}}{\partial\Delta\boldsymbol{\varepsilon}_{(n+1)}^{\mathrm{tr}}}$) is formularized following the procedure in [17,18].

$$\mathbf{D}^{ep} = \frac{\partial \Delta \tilde{\boldsymbol{\sigma}}_{(n+1)}}{\partial \Delta \boldsymbol{\varepsilon}_{(n+1)}^{\mathrm{tr}}} = \left[ (\mathbf{D}^e)^{-1} + \Delta \lambda \mathbf{P} + \frac{1}{2 \bar{\sigma}_y H \eta} \mathbf{M} \mathbf{n}_F \mathbf{n}_F^T - \frac{\Delta \lambda}{\eta^2} \mathbf{M} \mathbf{n}_F \mathbf{n}_F^T \mathbf{Z} \mathbf{P} \right]^{-1} \quad (24)$$

$$\text{where } \eta = \sqrt{\mathbf{n}_F^T \mathbf{Z} \mathbf{n}_F}, \quad H = \alpha\beta\tilde{p}_{(n+1)}^{\alpha-1}, \quad \mathbf{S}^e = (\mathbf{D}^e)^{-1}.$$

It is emphasized that the above operator is non-symmetric.

Therefore, a non-symmetric solver shall be invoked during defining the step in ABAQUS. It is worth mentioning that the quadratic formulation as presented above is very efficient to use in case of other plasticity functions or damage effect tensors by only changing the mapping matrices $\mathbf{P}$ and $\mathbf{Z}$ accordingly.
## 3 FEM IMPLEMENTATION

The detailed formulation of Ladevèze model presented in the preceding section has been implemented in ABAQUS using the user defined subroutines UMAT. The sequence of the implementation has been outlined in the flow chart referred in Fig.1.

<figure><img src="image_2.png" /><figcaption>**Fig.1**
UMAT subroutine implementation.</figcaption></figure>
# 4 RESULTS AND DISCUSSION

## 4.1 Elastic properties and model parameters

In the current work, material data reported in [8] were used to predict the results and compared with the experiments [20,21]. The elastic properties for Carbon Fiber Reinforced Polymers (CFRPs) and the strengths are presented in

Table 1., where $X_T$, $S_{12}$, $\varepsilon_{11}^f$ represent tensile failure strength in fiber direction, in-plane failure shear strength and failure strain in fiber direction respectively. The model specific parameters required for the Ladevèze model implementation are listed in Table 2., which were briefly discussed previously. Same lay-ups and dimensions of the composite coupons were used in the FEM simulation as were analyzed in [8] and are given in Table 3. A linear quadrilateral shell element S4R was used in all the simulations in the current work having element mesh size of 1 mm. On integration point was selected in the thickness direction based on the assumption that damage is same in out-of-plane direction in a single ply. Moreover, quarter symmetric geometry was modeled in ABAQUS to reduce the finite element domain and cost of computation.

<table><caption>**Table 1**
Elastic properties for the CFRP [8].</caption><thead><tr><th>E<sup>0</sup><sub>11</sub> (MPa)</th><th>E<sup>0</sup><sub>22</sub> (MPa)</th><th>G<sup>0</sup><sub>12</sub> (MPa)</th><th>ν<sub>12</sub></th><th>X<sub>T</sub> (MPa)</th><th>S<sub>12</sub> (MPa)</th><th>ε<sup>f</sup><sub>11</sub> (%)</th></tr></thead><tbody><tr><td>139000</td><td>10900</td><td>6000</td><td>0.32</td><td>2170</td><td>83</td><td>1.48</td></tr></tbody></table>

<table><caption>**Table 2**
Ladevèze model identification parameters[8].</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value (CFRP)</th></tr></thead><tbody><tr><td>Shear damage initiation threshold</td><td>Y<sub>12(0)</sub></td><td>0.048√MPa</td></tr><tr><td>Shear damage saturation level</td><td>Y<sub>12(c)</sub></td><td>3.10√MPa</td></tr><tr><td>Transvers damage initiation threshold</td><td>Y<sub>2(0)</sub></td><td>0.07√MPa</td></tr><tr><td>Transvers damage saturation level</td><td>Y<sub>2(c)</sub></td><td>2.75√MPa</td></tr><tr><td>Transvers brittle damage threshold</td><td>Y<sub>s</sub></td><td>0.565√MPa</td></tr><tr><td>Transverse-shear damage coupler</td><td>b</td><td>0.53</td></tr><tr><td>Transverse-shear plasticity coupler</td><td>a<sup>2</sup></td><td>0.54</td></tr><tr><td>Initial yield limit</td><td>σ̄<sub>0</sub></td><td>21.59 MPa</td></tr><tr><td>isotropic hardening law coefficient</td><td>β</td><td>558 MPa</td></tr><tr><td>isotropic hardening law exponent</td><td>α</td><td>0.38</td></tr></tbody></table>

<table><caption>**Table3**
Laminate configurations [8].</caption><thead><tr><th>Stacking sequence</th><th>L (mm)</th><th>w (mm)</th><th>t (mm)</th></tr></thead><tbody><tr><td>[0°]<sub>8</sub></td><td>250</td><td>15</td><td>1</td></tr><tr><td>[90°]<sub>16</sub></td><td>175</td><td>25</td><td>2</td></tr><tr><td>[±45°]<sub>28</sub></td><td>250</td><td>25</td><td>1</td></tr></tbody></table>

## 4.2 Model Prediction

Most of the UD FRPs with epoxy matrix behaves as linear-elastic with brittle fracture when the tensile load is applied along the fiber axis as the experiments showed in [8,22]. However, this behavior is also influenced by the fiber volume fraction, matrix material, fiber material and manufacturing processes reported in [14]. Thus excluding the plasticity in fiber direction, the failure criterion is simply either based on maximum normal stress theory or maximum strain theory which is not affected by the damage accumulation in the matrix. In the current model, the ultimate longitudinal failure strain is used proposed in [8,21].
$$ 
\left.
\begin{array}{l}
\varepsilon_{11} < \varepsilon_{11}^i \quad \text{then } d_1 = 0 \\
\\
\text{else if } \varepsilon_{11}^i < \varepsilon_{11} < \varepsilon_{11}^u \quad \text{then } d_1 = d_1^u \frac{\varepsilon_{11} - \varepsilon_{11}^i}{\varepsilon_{11}^u - \varepsilon_{11}^i} \\
\\
\text{else } d_1 = 1 - (1 - d_1^u) \frac{\varepsilon_{11}^u}{\varepsilon_{11}}
\end{array}
\right\} (25)
$$

In Eq.(25) $\varepsilon_{11}^i$ is the failure strain of the laminate in the fiber direction which is 1.48% in the present case. Any arbitrary value of $\varepsilon_{11}^u$ can be selected that must be greater than $\varepsilon_{11}^i$. The advantage of using this formulation for the fiber direction damage is that it avoids the numerical singularity problem and gives an asymptotic value of the damage rather than the abrupt loss of stiffness when the fiber strain reaches its failure strain. The predicted fiber-direction quasi-static behavior is depicted in Fig.2 (a).

<figure><img src="image_3.png" /><figcaption>**Fig.2** Numerical results comparison (a) Tension in fiber direction (b) Tension in Transverse direction.</figcaption></figure>

It can be noticed that this laminate behaves as linear elastic according to the experimental results taken from [8]. In addition, Fig.2 (b) demonstrates the transverse behavior to the fiber direction which is almost linear because plasticity is normally less dominated in this direction in CFRPs. Failure of FRPs in transverse tension loading is usually elastic and brittle.

In contrary to the fiber direction and transverse direction, a pronounced non-linearity is demonstrated in Fig. 3. All the CFRPs and glass fiber reinforced polymers (GFRPs) exhibit such non-linearity when the loading angle is near to 45° and is the consequence of the collaborated effect of the diffused damage and matrix plasticity. The accuracy of the model is evident as compared to the experiment. An elasto-plastic model based on the same plasticity potential has also been plotted in Fig. 3 where no damage coupling was considered in the UMAT. Evidently, the role of the coupled damage model with plasticity is considered indispensable to seek accurate modeling tools in order to exploit the composite to its limits in industrial applications.

<figure><img src="image_4.png" /><figcaption>**Fig.3** Shear stress-strain response.</figcaption></figure>
In order to gauge the robustness of the model, any arbitrary orientation laminates can be selected as the test cases. In [1] angle-ply laminates of [67.5°, 22.5°]₂s and [-12°, 78°]₂s were used as test cases whereas the mechanical behavior of angle-ply laminates having orientation of [±55°]₂s and [±80°]₂s were studied in [9] apart from the laminate [±67.5°]₂s which is also required for the determination of the transverse damage evolution law and both the plasticity and damage coupling parameters: *a*², *b*. In principle, the values of the aforementioned two coupling parameters shall be the same whether the [±67.5°]₂s is used or any other angle-ply laminate is used in the experiments where there prevail non-zero shear and transverse stresses, but the [±67.5°]₂s possesses some inherent edge over the others laminates. Both the transverse and shear stresses exist significantly in this laminate and thus produces a good coupling and at the same time gives less scatter in the data as emphasized in [1]. Despite the accuracy of the model, it requires unconventional testing for identification which is a drawback of this model as compared to the damage models based on the strength values and critical energy release rates for FRPs such as proposed in [16,18].

## 5 CONCLUSIONS

Keeping in view the importance and usefulness of the Ladevèze model in the composite community, a very comprehensive formulation of the model was presented. Particular attention was paid to the plasticity part of the model in quadratic form that is readily adaptable to use other plasticity potentials. The model was then implemented in the user material subroutines UMAT in ABAQUS with the derivation of the consistent elastoplastic operator. As based on the correlation of the prediction with the experimental results, it can be concluded that this model is comparatively an accurate and robust prediction tool. The limiting factor of the model is that when the base material (fiber, matrix) or the fiber volume fraction, and or the environmental conditions (temperature, humidity) are changed then the model identification parameters shall be determined experimentally for that new material system in new conditions. These model parameters need unconventional testing such as loading-unloading for the determination of the damage initiation and damage saturation which possesses the difficulty and somehow restrict this model to be introduced into the industrial environment. However, when the material is selected and once the parameters are extracted experimentally then this model is very appropriate to be used in the design offices for structure analysis of the FRPs.

## ACKNOWLEDGMENTS

The corresponding author, Israr Ud Din, highly acknowledges the financial support provided by Higher Education Commission (HEC), Pakistan in collaboration with Campus France for his doctoral studies based in France.

## REFERENCES

[1] Ladevèze PLD E., 1992, Damage modelling of the elementary ply for laminated composites, *Composites Science and Technology* **43**: 257-267.

[2] Liu PFZ J.Y., 2010, Recent developments on damage modeling and finite element analysis for composite laminates: A review, *Materials & Design* **31**(8): 3825-3834.

[3] Wang Y., Tong M., Zhu S., 2009, Three dimensional continuum damage mechanics model of progressive failure analysis in fibre-reinforced composite laminates, *Paper Presented at the Structures, Structural Dynamics, and Materials Conference*, California.

[4] Edgren F., Mattsson D., Asp L.E., Varna J., 2004, Formation of damage and its effects on non-crimp fabric reinforced composites loaded in tension, *Composites Science and Technology* **64**(5): 675-692.

[5] Herakovich C., Schroedter R., Gasser A., Guitard L., 2000, Damage evolution in [±45]s laminates with fiber rotation, *Composite Science and Technology* **60**: 2781-2789.

[6] Abu Al-Rub R.K.V., George Z., 2003, On the coupling of anisotropic damage and plasticity models for ductile materials, *International Journal of Solids and Structures* **40**(11): 2611-2643.

[7] Maire J., Chaboche J., 1997, A new formulation of continuum damage mechanics for composite materials, *Aeropace Science and Technology* **4**: 247-257.

[8] O'Higgins R.M., McCarthy C.T., McCarthy M.A., 2011, Identification of damage and plasticity parameters for continuum damage mechanics modelling of carbon and glass fibre-reinforced composite materials, *Strain* **47**(1): 105-115.
[9] Payan J., Hochard C., 2002, Damage modelling of laminated carbon/epoxy composites under static and fatigue loadings, *International Journal of Fatigue* **24**: 299-306.

[10] Greve L., Pickett A.K., 2006, Modelling damage and failure in carbon/epoxy non-crimp fabric composites including effects of fabric pre-shear, *Composites Part A: Applied Science and Manufacturing* **37**(11): 1983-2001.

[11] Boutaous A., Peseux B., Gornet L., Bélaidi A., 2006, A new modeling of plasticity coupled with the damage and identification for carbon fibre composite laminates, *Composite Structures* **74**(1): 1-9.

[12] Kachanov L.M., 1958, Time of the rupture process under creep condition, *Jzvestia Akad Nauk* **8**: 26-31.

[13] Chaboche J., 1988, Continuum damage mechanics: Part 1-General concepts, *Journal of Applied Mechanics* **55**: 59-64.

[14] Ribeiro M.L., Tita V., Vandepitte D., 2012, A new damage model for composite laminates, *Composite Structures* **94**(2): 635-642.

[15] Malgioglio F., Payan D., Magneville B., Farkas L., 2017, Material parameter identification challenge and procedure for intra-laminar damage prediction in unidirectional CFRP, *Paper Presented at the 6th ECCOMAS Thematic Conference on the Mechanical Response of Composites*.

[16] Chen J.F., Morozov E.V., Shankar K., 2012, A combined elastoplastic damage model for progressive failure analysis of composite materials and structures, *Composite Structures* **94**(12): 3478-3489.

[17] Neto E.A.S., Peric D., Owen D.R.J., 2008, *Computational Methods for Plasticity: Theory and Applications*, John Wiley and Sons, West Sussex, United Kingdom.

[18] Ud Din I., Hao P., Franz G., Panier S., 2018, Elastoplastic CDM model based on Puck's theory for the prediction of mechanical behavior of Fiber Reinforced Polymer (FRP) composites, *Composite Structures* **201**: 291-302.

[19] Simo J., Taylor R., 1985, Consistent tangent operators for rate-independent elastoplasticity, *Computer Methods in Applied Mechanics and Engineering* **48**: 101-118.

[20] O'Higgins R.M., McCarthy C.T., McCarthy M.A., 2012, Effects of shear-transverse coupling and plasticity in the formulation of an elementary Ply composites damage model, Part II: material characterization, *Strain* **48**(1): 59-67.

[21] O'Higgins R.M., McCarthy C.T., McCarthy M.A., 2012, Effects of shear-transverse coupling and plasticity in the formulation of an elementary Ply composites damage model, Part I: model formulation and validation, *Strain* **48**(1): 49-58.

[22] Truong T.C., Vettori M., Lomov S., Verpoest I., 2005, Carbon composites based on multi-axial multi-ply stitched preforms, *Applied Science and Manufacturing* **36**(9): 1207-1221.