# A viscoelastic-viscoplastic model for short-fibre reinforced polymers with complex fibre orientations

Mariem Nciri, Delphine Notta-Cuvier, Franck Lauro, Fahmi Chaari, Bassem Zouari, Yamen Maalej

► To cite this version:

Mariem Nciri, Delphine Notta-Cuvier, Franck Lauro, Fahmi Chaari, Bassem Zouari, et al.. A viscoelastic-viscoplastic model for short-fibre reinforced polymers with complex fibre orientations. DYMAT 2015 - 11th International Conference on the Mechanical and Physical Behaviour of Materials under Dynamic Loading, Sep 2015, Lugano, Switzerland. pp.04008, 10.1051/epjconf/20159404008 . hal-03433232

HAL Id: hal-03433232

https://uphf.hal.science/hal-03433232

Submitted on 29 Apr 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
# A viscoelastic-viscoplastic model for short-fibre reinforced polymers with complex fibre orientations

M. Nciri<fnref n="1" />,<fnref n="2" />,ª, D. Notta-Cuvier<fnref n="1" />, F. Lauro<fnref n="1" />, F. Chaari<fnref n="1" />, B. Zouari<fnref n="2" />, and Y. Maalej<fnref n="3" />

**Abstract.** This paper presents an innovative approach for the modelling of viscous behaviour of short-fibre reinforced composites (SFRC) with complex distributions of fibre orientations and for a wide range of strain rates. As an alternative to more complex homogenisation methods, the model is based on an additive decomposition of the state potential for the computation of composite's macroscopic behaviour. Thus, the composite material is seen as the assembly of a matrix medium and several linear elastic fibre media. The division of short fibres into several families means that complex distributions of orientation or random orientation can be easily modelled. The matrix behaviour is strain-rate sensitive, i.e. viscoelastic and/or viscoplastic. Viscoelastic constitutive laws are based on a generalised linear Maxwell model and the modelling of the viscoplasticity is based on an overstress approach. The model is tested for the case of a polypropylene reinforced with short-glass fibres with distributed orientations and subjected to uniaxial tensile tests, in different loading directions and under different strain rates. Results demonstrate the efficiency of the model over a wide range of strain rates.

## 1. Introduction

Thermoplastics reinforced with short fibres are more and more appealing for a wide-range of industrial applications, in particular in automotive industry in parts possibly subjected to sever loading conditions (e.g. crash...). As a consequence of complex reinforcement configurations in terms of orientation, geometrical and mechanical properties and complex polymeric matrix behaviour, SFRC behaviour is governed by numerous interdependent phenomena, such as plastic flow, ductile damage and strain-rate dependency of the matrix, fibre breakage, debonding or pull-out... Beyond the limits of linear elasticity, these phenomena can obviously not be modelled by a simple rule of mixture and more elaborated models, like the extended homogenization methods to components' non-linear behaviour have been developed. However, when all or some of phenomena listed above are interfering, the aforementioned methods become too complex to implement to allow to take all these relatively sophisticated features into account. To overcome these limitations, an original approach was recently proposed by Notta-Cuvier et al. [1,2], where complex fibre orientations can be easily dealt with thanks to the consideration of an assembly of the matrix medium with several fibre media, each of them being characterized by a volume fraction, $v_f^\alpha$, and a unit orientation vector. Stress transfer from the matrix to fibres is governed by a modified shear-lag model. Note that the model was then extended to strongly anisotropic ductile damage in the matrix material [3]. In this work, the

modelling is now extended to take the rate-dependency of the matrix behaviour into account.

## 2. Constitutive model

### 2.1. Stress state of matrix material

Rate dependency of the matrix behaviour is predicted with coupled viscoelastic-viscoplastic schemes. The viscoelastic part of the response is taken into account using the linear Wiechert model (i.e. Maxwell Generalized) consisting of a finite number of Maxwell elements arranged in parallel with a linear elastic Hooke element. This response is updated, when necessary by a viscoplastic scheme. The hereditary integral over the time, t, of the viscoelastic stress response is formulated as:

$$ \sigma(t) = \int_{-\infty}^{t} R^{ve}(t-\Gamma) : \frac{d\varepsilon(\zeta)}{d\zeta} d\zeta \quad (1) $$

where $\sigma$ and $\varepsilon$ are the matrix stress and strain tensors, respectively. $R^{ve}$ is the fourth order relaxation tensor expressed by:

$$ R^{ve}(t) = L_{\infty}^{ve} + \sum_{i=1}^{N} L_{i}^{ve} \exp \left( -\frac{t}{\tau_{i}} \right) \quad (2) $$

where N is the number of Maxwell elements. $L_{\infty}^{ve}$ is the fourth order long-term elasticity tensor of the Hooke's element defined by:

$$ L_{\infty}^{ve} = 2G_{\infty}I_{d} + K_{\infty}I \otimes I \quad (3) $$

<footnote n="1">¹University of Valenciennes (UVHC), LAMIH UMR CNRS 8201, 59313 Valenciennes, France</footnote>

<footnote n="2">²National Engineering School of Sfax (ENIS), L2MP, BP. W3038, Sfax, Tunisia</footnote>

<footnote n="3">³ University of Tunis El Manar, ENIT, MAI (LR11ES19), 1002 Tunis, Tunisia</footnote>

<footnote>ª Corresponding author: mariem.nciri@etu.univ-valenciennes.fr</footnote>

<footnote>This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.</footnote>
$\tau_i$ and $L_i^{ve}$ are the relaxation time and the fourth order elasticity tensor of the ith Maxwell element, respectively, with:

$$L_i^{ve} = 2G_i I_d + K_i I \otimes I \quad (4)$$

$I_d$ is the deviatoric projection tensor such as:

$$I_d = I_s - \frac{1}{3} I \otimes I. \quad (5)$$

The matrix viscoplasticity is taken into account following Perzyna's viscoplastic model [4]. The pressure dependent and nonsymmetric (visco) plastic flow, i.e. sensitivity of the polymeric matrix behaviour to the nature of loading (e.g., tension, compression..), are modelled by Raghava's yield surface [5]. Plastic flow occurs as soon as the first invariant of the stress tensor, $I_1(\sigma)$, and the second invariant of the deviatoric stress tensor, $I_2(\sigma)$, reach a critical combination. Expression of the yield surface, $f$, is given by:

$$f(\sigma, R) = \frac{(\eta - 1) I_1(\sigma) + \sqrt{(\eta - 1)^2 I_1^2(\sigma) + 12\eta I_2(\sigma)}}{2\eta} - \sigma_t - R(\kappa). \quad (6)$$

with

$$I_1(\sigma) = tr(\sigma) \quad (7a)$$

and

$$I_2(\sigma) = \frac{1}{2} S : S \quad (7b)$$

where $S$ is the deviatoric stress tensor. The hydrostatic pressure dependency parameter, $\eta$, is defined using the ratio between the quasi-static initial yield stress in compression and tension, $\sigma_c$ and $\sigma_t$, respectively, and:

$$\eta = \frac{\sigma_c}{\sigma_t}. \quad (8)$$

In the expression of the yield surface, $R(\kappa)$, is a nonlinear isotropic hardening function, defined by Balieu et al. [6]:

$$R(\kappa) = Q_1 \kappa \exp(-b_1 \kappa) + Q_2(1 - \exp(-b_2 \kappa)) + b_3 \kappa^3 + b_4 \kappa^2 + b_5 \kappa \quad (9)$$

where $Q_1$, $Q_2$, $b_1$, $b_2$, $b_3$, $b_4$ and $b_5$ are material parameters. The hardening variable $\kappa$ is the equivalent viscoplastic strain defined by:

$$\kappa = \sqrt{\frac{3}{2} \varepsilon^{vp} : \varepsilon^{vp}} \quad (10)$$

with $\varepsilon^{vp}$ the viscoplastic strain tensor. Since the viscoplastic deformation of polymer materials is not an isochoric phenomenon, a non-associated viscoplastic flow rule is used. The asymmetry of this volume change (i.e., not similar for positive and negative hydrostatic pressures) is modelled by a hyperbolic viscoplastic dissipation potential [6], defined by:

$$F(\sigma) = \sqrt{3I_2(\sigma) + \frac{1}{3}(a^+ < I_1 >^2 + a^- < -I_1 >^2)} \quad (11)$$

where $a^+$ and $a^-$ are volume variation parameters for positive and negative hydrostatic pressures, respectively. The viscoplastic strain rate tensor can therefore evolve differently for positive and negative hydrostatic pressures. The symbol $\langle \cdot \rangle$ is the Macauley braket, that is, for any scalar $x$, given by $\frac{(x+|x|)}{2}$. According to this non-associated flow rule the viscoplastic strain rate tensor is given by:

$$\dot{\varepsilon}^{vp} = \dot{\lambda}n. \quad (12)$$

Where $n$ is the viscoplastic flow direction tensor expressed by:

$$n = \frac{\partial F}{\partial \sigma}. \quad (13)$$

Consequently, with the given dissipation potential, the viscoplastic strain rate tensor becomes:

$$\dot{\varepsilon}^{vp} = \dot{\lambda} \frac{\frac{3}{2} S + \frac{1}{9} (a^+ < I_1 > + a^- < -I_1 >) I}{\sqrt{3I_2(\sigma) + \frac{1}{27} (a^+ < I_1 >^2 + a^- < -I_1 >^2)}} \quad (14)$$

Matrix material's strain rate dependency is introduced by a viscoplastic formulation. The static yield surface, $f$, is extended according to an overstress based viscoplastic theory and a dynamic yield surface is then defined, as follows [4]:

$$F^{vp}(\sigma, R, \dot{\kappa}) = \frac{(\eta - 1) I_1(\sigma) + \sqrt{(\eta - 1)^2 I_1^2(\sigma) + 12\eta I_2(\sigma)}}{2\eta} - (\sigma_t + R(\kappa)) - \dot{\sigma}^{vp} \equiv 0. \quad (15)$$

With $\dot{\sigma}^{vp}$ the overstress: As postulated in Perzyna's model, the viscoplastic multiplier $\dot{\lambda}$ is expressed in terms of overstress as follows:

$$\dot{\lambda} = \frac{\dot{\kappa}_0}{\sqrt{\frac{2}{3}n : n}} \left( \frac{\dot{\sigma}^{vp}}{2\eta(\sigma_t + R(\kappa))} \right)^{\frac{1}{n}} \quad (16)$$

and equivalently:

$$\dot{\sigma}^{vp} = (\sigma_t + R(\kappa)) \left( \frac{\dot{\kappa}}{\dot{\kappa}_0} \right)^n \quad (17)$$

where $n$ and $\dot{\kappa}_0$ are the strain rate sensitivity and viscosity parameters, respectively, and $\dot{\kappa}$ is the equivalent viscoplastic strain rate, defined by:

$$\dot{\kappa} = \sqrt{\frac{2}{3} \dot{\varepsilon}^{vp} : \dot{\varepsilon}^{vp}} = \dot{\lambda} \sqrt{\frac{2}{3} n : n} \quad (18)$$

The expression of the dynamic yield surface can therefore be rewritten as follows:

$$F^{vp}(\sigma, R, \dot{\kappa}) = \frac{(\eta - 1) I_1(\sigma) + \sqrt{(\eta - 1)^2 I_1^2(\sigma) + 12\eta I_2(\sigma)}}{2\eta} - (\sigma_t + R(\kappa)) \left(1 + \frac{\dot{\kappa}}{\dot{\kappa}_0}\right)^n \equiv 0. \quad (19)$$
The standard Kuhn-Tucker loading/unloading conditions (i.e., $F^{vp} \le 0$, $\dot{\lambda} \ge 0$, $\dot{\lambda}F^{vp} = 0$) are applied to the dynamic yield surface. It can be noted that both static and dynamic yield surfaces are updated all along the implementation. As a consequence, during unloading overstress does not systematically vanish and stress state can remain in the viscoplastic domain (i.e. stress state above the static yield surface, on updated dynamic surface). The rate form of the viscoplastic multiplier is obtained by substituting its expression into the viscous yield surface and is as follows:

$$ \dot{\lambda} = \begin{cases} 0, & \text{if } f < 0, \\ \frac{\kappa_0}{\sqrt{\frac{2}{3}n:n}} \left( \frac{(\eta-1)I_1(\sigma) + \sqrt{(\eta-1)^2 I_1^2(\sigma) + 12\eta I_2(\sigma)}}{2\eta(\sigma_1 + R(\kappa))} \right)^{\frac{1}{n}}, & \text{if } f \ge 0. \end{cases} \quad (20) $$

## 2.2. Stress state of the composite material

Once matrix behaviour is solved, contributions of fibres media to SFRC behaviour are introduced according to the original method developed by Notta-Cuvier et al. [1]. Each fibre family is subjected to the projection of the total deformation gradient tensor along its orientation and behaves in a one-dimensional linear elastic way. The average axial stress in the fibres is computed based on Bowyer-Bader theory [7], where a progressive degradation of the fibre-matrix interface is taken into account. As previously stated, the stress state of the composite material is computed as a combination of the different constituent contributions, starting from a state potential, the Helmholtz free energy (Eq. (21)), that is divided into as many terms as there are media constituting the composite.

$$ \rho\phi = \nu_M\rho_M\phi_M + \sum_{\alpha=1}^{n_{fam}} \nu_F^\alpha \rho_F^\alpha \phi_F^\alpha \quad (21) $$

where $\rho$, $\rho_M$ and $\rho_F^\alpha$ are respectively the densities of the composite material, the matrix material and the fibres family $\alpha$. Using a thermodynamic formulation a particular solution of the Clausius-Duhem inequality can be expressed by Eq. (21) thus giving the expression of composite macroscopic stress tensor.

$$ \sigma = \nu_M\sigma_M + \sum_{\alpha=1}^{n_{fam}} \nu_F^\alpha A^\alpha \sigma_F^\alpha A^\alpha. \quad (22) $$

## 3. Example of application

The model is implemented as a VUMAT (Abaqus Explicit) and is tested for a fictitious injection moulded SFRC, namely a polypropylene reinforced with 10 vol.% glass fibres. Fibre configuration is extracted from work by Kammoun et al. [8] ($E_F = 76$ GPa, distributed orientations with preferential orientation in the injection flow direction). In the present implementation, the fibre distribution is modelled using 10 fibre media with angles of orientation varying from 0 to 90°, by step of 10°, with respect to the injection flow direction and a suitable

<figure><img src="image_2.png" /><figcaption>Figure 1. Distribution of orientations [8].</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. $\theta = 0^\circ$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. $\theta = 30^\circ$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. $\theta = 45^\circ$.</figcaption></figure>

volume fraction is associated to each fibre media, as shown in Fig. 1. The matrix viscoelastic-viscoplastic material parameters are chosen in accordance with the model proposed by Balieu et al. [6], for a 20% mineral filled polypropylene co-polymer.

In order to assess the model accuracy, uniaxial tensile tests are performed on one element for different strain rates and loading directions (angle of loading, $\theta$, successively equals to 0°, 30°, 45° and 90° with respect to injection flow direction) as shown in Figs. 2 to 5.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** $\theta = 90^{\circ}$.</figcaption></figure>

Responses' dependency on the strain rate is highlighted in the linear and non-linear phases according to the hereinabove presented implementation. This evolution is obviously dependent on the loading direction regarding the principal fibre direction (i.e. the injection flow direction). It can be seen that, as the angle of loading direction with respect to fibres' one increases, the strain rate dependency of the whole material's modelled response is dimmed down, which is in accordance with the fact that the rate independent linear elastic response of the fibres is of greater influence in its principal direction, while at increased angles of orientation the response is rather dominated by that of the matrix.

## 4. Conclusions

A behaviour model for strain-rate dependent short fibre reinforced composites is presented, based on an original approach that aims to be an efficient alternative to more complex procedures of homogenisation. The matrix has viscoelastic-viscoplastic and pressure dependent behaviour and complex reinforcement configurations are considered. Current development concerns a physical modelling of fibre/matrix interface

degradation during loading, with possible strain-rate effects. The determination of the material constants is the subject of an already started static and dynamic experimental investigation.

This research is conducted through collaboration between the University of Valenciennes and the National Engineering School of Sfax. This collaboration is jointly financed in the frame of the Unique CMCU programme. The present research work has been supported by International Campus on Safety and Intermodality in Transportation, the Région Nord Pas de Calais, the European Community, the Délégation Régionale à la Recherche et à la Technologie, the Ministère de l'Enseignement Supérieur et de la Recherche, and the Centre National de la Recherche Scientifique: the authors gratefully acknowledge the support of these institutions.

## References

[1] D. Notta-Cuvier, F. Lauro, B. Bennani, R. Balieu. International Journal of Solids and Structures **50**, 2857-2871 (2013)

[2] D. Notta-Cuvier, F. Lauro, B. Bennani. Composites: Part A **62**, 60-66 (2014)

[3] D. Notta-Cuvier, F. Lauro, B. Bennani, R. Balieu. Mechanics of materials **68**, 193-206 (2014)

[4] P. Perzyna. Advances in Applied Mechanics **9**, 243-377 (1966)

[5] R. Raghava, R.M. Caddell, G.S.Y. Yeh. Journal of Materials Sciences **8**, 225-232 (1973)

[6] R. Balieu, F. Lauro, B. Bennani, R. Delille, T. Matsumoto, E. Mottola. International Journal of Plasticity **51**, 241-271 (2013)

[7] W.H. Bowyer, M.G. Bader. Journal of Materials Science **7**, 1315-1321 (1972)

[8] S. Kammoun, I. Doghri, L. Adam, G. Robert, L. Delannay. Composites Part A **42**, 2892-2902 (2011)