# Anisotropic visco-elastoplastic modeling of quasi-unidirectional flax riber reinforced epoxy behavior: an investigation on low-velocity impact response

Marwa Abida, Jamel Mars, Florian Gehring, Alexandre Vivet, Fakhreddine Dammak

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

## To cite this version:

Marwa Abida, Jamel Mars, Florian Gehring, Alexandre Vivet, Fakhreddine Dammak. Anisotropic visco-elastoplastic modeling of quasi-unidirectional flax riber reinforced epoxy behavior: an investigation on low-velocity impact response. Journal of Renewable Materials, 2018, 6 (5), pp.464-476. 10.32604/JRM.2018.01897. hal-02339891

HAL Id: hal-02339891

https://hal.science/hal-02339891

Submitted on 28 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Anisotropic Visco-Elastoplastic Modeling of Quasi-Unidirectional Flax Fiber Reinforced Epoxy Behavior: An Investigation on Low-Velocity Impact Response

Marwa Abida<fnref n="1" />,<fnref n="2" />, Jamel Mars<fnref n="3" />,*, Florian Gehring<fnref n="1" />, Alexandre Vivet<fnref n="1" /> and Fakhreddine Dammak<fnref n="3" />

**ABSTRACT:** Based on experimental test results, flax fiber reinforced polymer composites are characterized by nonlinear visco-elastoplastic behavior. The aim of this work is to model the quasi-unidirectional flax fiber reinforced composite behavior through a three dimensional formulation with orthotropic elasticity and orthotropic plasticity using Hill criterion. The isotropic hardening and Johnson-Cook parameters are identified from unidirectional tensile tests at different strain rates. The adjustment of Hill's yield criterion is developed based on yield stresses obtained in tensile tests at different directions. The numerical integration of the constitutive equations is implemented in a user-defined material, UMAT subroutines for the commercial finite element code ABAQUS. Once model parameters are identified using tensile tests, the model needs to be validated by confronting it with other experimental results. That is why experimental and numerical three-point bending tests are carried out in order to validate the proposed model with tests that have not served for the identification. Finally, a numerical parametric study on low velocity impact of a flax/epoxy composite circular plate is investigated.

**KEYWORDS:** Visco-elastoplasticity, anisotropic behavior, impact simulation, flax fiber, composites

## 1 INTRODUCTION

Conventional composites are constituted of synthetic fibers, such as glass fibers or carbon fibers, reinforced with oil-based matrices like epoxy or polypropylene. In recent years, new environmental protocols and policies have encouraged the return of natural fibers in the forefront of the scene [1]. As they originate from renewable resources, plant fibers represent an interesting alternative to synthetic fibers [2].

Due to their low density, their reasonable cost coupled with their high specific mechanical properties, and their availability especially at the European scale, flax fibers can be considered as the most interesting vegetal fibers [3]. Their performance has been investigated for several years. Flax fibers are traditionally used in the textile industry. Recently, they have been employed in several other domains such as transport, leisure and sport for instance. Flax fibers can also be combined with mineral fibers to optimize the properties of the assembly [4].

In order to use flax fiber reinforced polymer (FFRP) composites in industrial and technical applications, a rough comprehension and a reliable prediction of

their behavior are needed. However, there is a dearth of information in this domain. A prediction of the non-linear stress/strain behavior of unidirectional and quasi-unidirectional flax fiber reinforced composites was funded via the development of a semi-empirical tensor linear model [5]. Rubio-López et al. [6] established a unidirectional model to predict the viscoplastic behavior of composites reinforced with natural fibers, which are flax, cotton and jute. A unidirectional constitutive model considering the viscoelastic behavior of flax fiber reinforced epoxy composites was proposed by Poilane et al. [7]. A non-linear phenomenological model based on rheological elements with eight parameters was adopted.

Three-dimensional modeling of the quasi-unidirectional flax/epoxy composites behavior is, according to the author's knowledge, poorly raised in the literature. The flax epoxy composites behavior is anisotropic in the elastic region or outside the elastic range. This anisotropy has been perceived experimentally via tensile tests at different fiber orientations [8]. An appropriate material model capable of being used in a reliable design must encompass the complex flax/epoxy composites behavior.

Nowadays, natural fiber reinforced composites are more and more used in automotive applications [9]. These FFRP composites are widely used in many

<footnote n="1">¹Normandie Univ, ENSICAEN, UNICAEN, CEA, CNRS, CIMAP, 14000 Caen, France</footnote>

<footnote n="2">²Mechanical Modeling and Manufacturing Laboratory (LA2MP), National School of Engineers of Sfax, University of Sfax, B.P. 1173-3038, Sfax, Tunisia</footnote>

<footnote n="3">³Laboratory of Electromechanical Systems (LASEM), National Engineering School of Sfax, University of Sfax, B.P. 1173-3038, Sfax, Tunisia</footnote>

<footnote>*Corresponding author: jamelmars@yahoo.fr</footnote>

<footnote>DOI: 10.32604/JRM.2018.01897</footnote>
# Anisotropic Visco-Elastoplastic Modeling of Quasi-Unidirectional Flax Fiber Reinforced Epoxy Behavior: An Investigation on Low-Velocity Impact Response

Marwa Abida<fnref n="1" />,<fnref n="2" />, Jamel Mars<fnref n="3" />,*, Florian Gehring<fnref n="1" />, Alexandre Vivet<fnref n="1" /> and Fakhreddine Dammak<fnref n="3" />

**ABSTRACT:** Based on experimental test results, flax fiber reinforced polymer composites are characterized by nonlinear visco-elastoplastic behavior. The aim of this work is to model the quasi-unidirectional flax fiber reinforced composite behavior through a three dimensional formulation with orthotropic elasticity and orthotropic plasticity using Hill criterion. The isotropic hardening and Johnson-Cook parameters are identified from unidirectional tensile tests at different strain rates. The adjustment of Hill's yield criterion is developed based on yield stresses obtained in tensile tests at different directions. The numerical integration of the constitutive equations is implemented in a user-defined material, UMAT subroutines for the commercial finite element code ABAQUS. Once model parameters are identified using tensile tests, the model needs to be validated by confronting it with other experimental results. That is why experimental and numerical three-point bending tests are carried out in order to validate the proposed model with tests that have not served for the identification. Finally, a numerical parametric study on low velocity impact of a flax/epoxy composite circular plate is investigated.

**KEYWORDS:** Visco-elastoplasticity, anisotropic behavior, impact simulation, flax fiber, composites

## 1 INTRODUCTION

Conventional composites are constituted of synthetic fibers, such as glass fibers or carbon fibers, reinforced with oil-based matrices like epoxy or polypropylene. In recent years, new environmental protocols and policies have encouraged the return of natural fibers in the forefront of the scene [1]. As they originate from renewable resources, plant fibers represent an interesting alternative to synthetic fibers [2].

Due to their low density, their reasonable cost coupled with their high specific mechanical properties, and their availability especially at the European scale, flax fibers can be considered as the most interesting vegetal fibers [3]. Their performance has been investigated for several years. Flax fibers are traditionally used in the textile industry. Recently, they have been employed in several other domains such as transport, leisure and sport for instance. Flax fibers can also be combined with mineral fibers to optimize the properties of the assembly [4].

In order to use flax fiber reinforced polymer (FFRP) composites in industrial and technical applications, a rough comprehension and a reliable prediction of

their behavior are needed. However, there is a dearth of information in this domain. A prediction of the nonlinear stress/strain behavior of unidirectional and quasi-unidirectional flax fiber reinforced composites was funded via the development of a semi-empirical tensor linear model [5]. Rubio-López et al. [6] established a unidirectional model to predict the viscoplastic behavior of composites reinforced with natural fibers, which are flax, cotton and jute. A unidirectional constitutive model considering the viscoelastic behavior of flax fiber reinforced epoxy composites was proposed by Poilane et al. [7]. A nonlinear phenomenological model based on rheological elements with eight parameters was adopted.

Three-dimensional modeling of the quasi-unidirectional flax/epoxy composites behavior is, according to the author's knowledge, poorly raised in the literature. The flax epoxy composites behavior is anisotropic in the elastic region or outside the elastic range. This anisotropy has been perceived experimentally via tensile tests at different fiber orientations [8]. An appropriate material model capable of being used in a reliable design must encompass the complex flax/epoxy composites behavior.

Nowadays, natural fiber reinforced composites are more and more used in automotive applications [9]. These FFRP composites are widely used in many

<footnote n="1">¹Normandie Univ, ENSICAEN, UNICAEN, CEA, CNRS, CIMAP, 14000 Caen, France</footnote>

<footnote n="2">²Mechanical Modeling and Manufacturing Laboratory (LA2MP), National School of Engineers of Sfax, University of Sfax, B.P. 1173-3038, Sfax, Tunisia</footnote>

<footnote n="3">³Laboratory of Electromechanical Systems (LASEM), National Engineering School of Sfax, University of Sfax, B.P. 1173-3038, Sfax, Tunisia</footnote>

<footnote>*Corresponding author: jamelmars@yahoo.fr</footnote>

<footnote>DOI: 10.32604/JRM.2018.01897</footnote>
engineering applications such as automobile and aerospace. During their manufacture, maintenance and service life, these FFRP composites may suffer different types of damage, of which low velocity impact is considered as one of the most important and dangerous, because it can induce internal damage leading to a significant strength reduction of the structure. However, these properties are difficult to quantify. Impact testing of materials is conducted to determine the amount of energy that can be absorbed by a material when a force is applied suddenly. Usually, impact testing is performed using Charpy and Izod test machines to determine ductile-brittle transitions in metals [10]. Charpy and Izod tests assume the pre-existence of a notch which is not suitable for testing composite materials [11]. Due to the absence of a notch, the rupture in composites cannot be controlled and the measured failure energy depends on several phenomena such as delamination, fiber failure, matrix failure and fiber matrix debonding in uncontrolled proportions. The drop-weight test is an alternative method, which consists in dropping a known mass from a given height onto a flat sample. As it expresses better what a structural component would experience in its service life, drop weight test is considered as a more realistic test which has proven to be a valuable source of information about the impact behavior of woven composites [12], tape laminates [13] and sandwich structures [14].

Low velocity impact behavior of composite laminates has been investigated intensively [15-19], among authors. However, limited works can be found on low velocity impact behavior of plant fiber reinforced composites. Recently, some works were conducted on low velocity impact behavior of hybrid composites, which consist of a combination of synthetic and plant fibers impregnated with a polymer matrix [20, 21]. Recent studies have begun to explore the impact behavior of hemp fiber reinforced composites [22-24]. A finite element numerical model was established by Rubio-López et al. [6] in order to predict the low velocity impact behavior of all-cellulose composite plates made from Cordenka woven plies. Liang et al., [25] conducted an experimental study of flax/epoxy impact behavior of composite plates subjected to low-velocity transverse impact.

In the present work, a three-dimensional material modeling which can be implemented in current finite element codes is presented. The model is able to describe the complex behavior of flax/epoxy composites. This model takes into account the orthotropic elastic behavior, the anisotropic plasticity

using Hill criterion (1948) [26], the isotropic hardening and the viscoplastic behavior. Tensile results at different strain rates conducted by Poilane et al. [7] are used to determine the viscoplastic behavior of the flax/epoxy composites. Elastic mechanical properties of the same composite obtained by Cherif et al. [8] are deployed to establish the three-dimensional anisotropic elasticity. Tensile tests at different fiber orientations [8] are used to extract Hill anisotropic coefficients. Three-point bending tests are carried out experimentally and numerically in order to validate the model with tests that have not served for the identification. The model is then used to simulate low velocity impact behavior of flax/epoxy composites. The present proposed three-dimensional formulation is the first attempt to model the quasi-unidirectional flax epoxy composite behavior with orthotropic elasticity and plasticity. The remainder of this paper is as follows. The constitutive equations development and the numerical integration are given in Section 2. In Section 3, model identification and validation of flax/epoxy composite is established. Experimental tests on quasi-unidirectional flax/epoxy composite are explored in order to establish a 3D viscoplastic with anisotropic elasticity model. Section 4 consists in modeling the impact behavior of epoxy/flax circular plate. A parametric study is accomplished based on the impactor mass, initial velocities, impactor geometry and the circular plate thickness.

## 2 CONSTITUTIVE EQUATIONS AND NUMERICAL IMPLEMENTATION

### 2.1 Three-Dimensional Rate-Independent Orthotropic Elasto-Plasticity

A three-dimensional rate-independent model with isotropic hardening, which takes into account the anisotropic behavior is established. The total strain tensor $\varepsilon$ is composed of an elastic strain tensor $\varepsilon^e$ (instantaneous reversible strain) and a plastic (time-independent irreversible) strain tensor $\varepsilon^p$.

$$ \varepsilon = \varepsilon^e + \varepsilon^p \quad (1) $$

The stress tensor $\sigma$ is related to the elastic strain by the elastic stiffness tensor $D$ as follows

$$ \sigma = D : \varepsilon^e \quad (2) $$

In the case of orthotropic elastic materials, the inverse of the elasticity matrix $D$, referred to as the compliance matrix $S$, is expressed as follows
engineering applications such as automobile and aerospace. During their manufacture, maintenance and service life, these FFRP composites may suffer different types of damage, of which low velocity impact is considered as one of the most important and dangerous, because it can induce internal damage leading to a significant strength reduction of the structure. However, these properties are difficult to quantify. Impact testing of materials is conducted to determine the amount of energy that can be absorbed by a material when a force is applied suddenly. Usually, impact testing is performed using Charpy and Izod test machines to determine ductile-brittle transitions in metals [10]. Charpy and Izod tests assume the pre-existence of a notch which is not suitable for testing composite materials [11]. Due to the absence of a notch, the rupture in composites cannot be controlled and the measured failure energy depends on several phenomena such as delamination, fiber failure, matrix failure and fiber matrix debonding in uncontrolled proportions. The drop-weight test is an alternative method, which consists in dropping a known mass from a given height onto a flat sample. As it expresses better what a structural component would experience in its service life, drop weight test is considered as a more realistic test which has proven to be a valuable source of information about the impact behavior of woven composites [12], tape laminates [13] and sandwich structures [14].

Low velocity impact behavior of composite laminates has been investigated intensively [15-19], among authors. However, limited works can be found on low velocity impact behavior of plant fiber reinforced composites. Recently, some works were conducted on low velocity impact behavior of hybrid composites, which consist of a combination of synthetic and plant fibers impregnated with a polymer matrix [20, 21]. Recent studies have begun to explore the impact behavior of hemp fiber reinforced composites [22-24]. A finite element numerical model was established by Rubio-López et al. [6] in order to predict the low velocity impact behavior of all-cellulose composite plates made from Cordenka woven plies. Liang et al., [25] conducted an experimental study of flax/epoxy impact behavior of composite plates subjected to low-velocity transverse impact.

In the present work, a three-dimensional material modeling which can be implemented in current finite element codes is presented. The model is able to describe the complex behavior of flax/epoxy composites. This model takes into account the orthotropic elastic behavior, the anisotropic plasticity

using Hill criterion (1948) [26], the isotropic hardening and the viscoplastic behavior. Tensile results at different strain rates conducted by Poilane et al. [7] are used to determine the viscoplastic behavior of the flax/epoxy composites. Elastic mechanical properties of the same composite obtained by Cherif et al. [8] are deployed to establish the three-dimensional anisotropic elasticity. Tensile tests at different fiber orientations [8] are used to extract Hill anisotropic coefficients. Three-point bending tests are carried out experimentally and numerically in order to validate the model with tests that have not served for the identification. The model is then used to simulate low velocity impact behavior of flax/epoxy composites. The present proposed three-dimensional formulation is the first attempt to model the quasi-unidirectional flax epoxy composite behavior with orthotropic elasticity and plasticity. The remainder of this paper is as follows. The constitutive equations development and the numerical integration are given in Section 2. In Section 3, model identification and validation of flax/epoxy composite is established. Experimental tests on quasi-unidirectional flax/epoxy composite are explored in order to establish a 3D viscoplastic with anisotropic elasticity model. Section 4 consists in modeling the impact behavior of epoxy/flax circular plate. A parametric study is accomplished based on the impactor mass, initial velocities, impactor geometry and the circular plate thickness.

## 2 CONSTITUTIVE EQUATIONS AND NUMERICAL IMPLEMENTATION

### 2.1 Three-Dimensional Rate-Independent Orthotropic Elasto-Plasticity

A three-dimensional rate-independent model with isotropic hardening, which takes into account the anisotropic behavior is established. The total strain tensor $\epsilon$ is composed of an elastic strain tensor $\epsilon^e$ (instantaneous reversible strain) and a plastic (time-independent irreversible) strain tensor $\epsilon^p$.

$$ \epsilon = \epsilon^e + \epsilon^p \quad (1) $$

The stress tensor $\sigma$ is related to the elastic strain by the elastic stiffness tensor $D$ as follows

$$ \sigma = D : \epsilon^e \quad (2) $$

In the case of orthotropic elastic materials, the inverse of the elasticity matrix $D$, referred to as the compliance matrix $S$, is expressed as follows
$$ S = \begin{bmatrix} D = S^{-1} \\ \begin{bmatrix} \frac{1}{E_1} & -\frac{\nu_{12}}{E_1} & -\frac{\nu_{13}}{E_1} & 0 & 0 & 0 \\ \frac{1}{E_2} & -\frac{\nu_{23}}{E_2} & 0 & 0 & 0 & 0 \\ \frac{1}{E_3} & 0 & 0 & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{\frac{1}{G_{12}} & 0 & 0} \\ \multicolumn{6}{c}{\frac{1}{G_{13}} & 0} \\ \multicolumn{6}{c}{\frac{1}{G_{23}}} \end{bmatrix} \end{bmatrix}, \quad (3) $$

where $(E_1, E_2, E_3)$ are longitudinal and transverse Young's moduli, $(G_{12}, G_{13}, G_{23})$ are shear moduli and $(\nu_{12}, \nu_{13}, \nu_{23})$ are Poisson ratios. With the plane stress assumption, the elasticity matrix $D$ became

$$ D = \begin{bmatrix} H_{11} & H_{12} & 0 & 0 & 0 \\ H_{22} & 0 & 0 & 0 & 0 \\ 0 & G_{12} & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{G_{13}} \\ \multicolumn{6}{c}{G_{23}} \end{bmatrix} \quad (4) $$

where $G_{12} = G_{13}$. The components of $D$ can be written in terms of Engineering constants as follows:

$$ \left\{ \begin{array}{l} H_{11} = \frac{E_1}{a}, H_{22} = \frac{E_2}{a}, H_{12} = \frac{\nu_{12}E_2}{a} = \frac{\nu_{21}E_1}{a} \\ a = 1 - \nu_{12}\nu_{21} = 1 - \frac{\nu_{12}^2 E_2}{E_1} \end{array} \right. \quad (5) $$

Concerning the yield criterion, it is chosen as

$$ f = \varphi(\boldsymbol{\sigma}) - \sigma_p, \quad \sigma_p = \sigma_Y + R \quad (6) $$

where $\sigma_Y$ is the initial yield stress, $\varphi(\boldsymbol{\sigma})$ is the equivalent stress and $R$ is the drag stress in isotropic hardening. In order to express the anisotropic plasticity, the equivalent stress $\varphi(\boldsymbol{\sigma})$ can be described by using the Hill criterion [26]

$$ \varphi(\boldsymbol{\sigma}) = \|\boldsymbol{\sigma}\|_{H} = \sqrt{\boldsymbol{\sigma}^{T} \boldsymbol{H} \boldsymbol{\sigma}} \quad (7) $$

where $H$ is a fourth-order tensor given in a condensed form as [27-29].

$$ H = \begin{bmatrix} H+G & -H & -G & 0 & 0 & 0 \\ H+F & -F & 0 & 0 & 0 & 0 \\ F+G & 0 & 0 & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{2N} \\ \multicolumn{6}{c}{2M} \\ \multicolumn{6}{c}{2L} \end{bmatrix} \quad (8) $$

$H, G, F, N, M$ and $L$ are Hill material constants to be identified. The evolution equations for the present problem are given as follows

$$ \dot{\boldsymbol{\epsilon}}^p = \dot{\gamma} \frac{\partial f}{\partial \boldsymbol{\sigma}} = \dot{\gamma} \boldsymbol{n}, \quad \boldsymbol{n} = \frac{1}{\varphi(\boldsymbol{\sigma})} \boldsymbol{H} \boldsymbol{\sigma} \quad (9) $$

$$ \dot{r} = -\dot{\gamma} \frac{\partial f}{\partial R} = \dot{\gamma} \quad (10) $$

where $\dot{\gamma}$ is the plastic multiplier which is consistent with the loading/unloading condition

$$ \dot{\gamma} \ge 0, \quad f \le 0, \quad \dot{\gamma} f = 0 \quad (11) $$

## 2.2 Numerical Integration

A fully implicit integration algorithm is used to solve the constitutive equations. The main advantage of the present local fully implicit integration algorithm is its strain driven characteristic and unconditional stability with a low number of equations to be solved: only one equation. The unknown of this equation is the plastic multiplier $\Delta\gamma$

$$ f(\Delta\gamma) = \varphi(\boldsymbol{\sigma}) - \sigma_p = 0, \quad \varphi = \left[ \boldsymbol{\sigma}^{\mathrm{tr}T} \boldsymbol{I}_c^{-T} \boldsymbol{H} \boldsymbol{I}_c^{-1} \boldsymbol{\sigma}^{\mathrm{tr}} \right]^{1/2} \quad (12) $$

$$ \boldsymbol{\sigma}^{\mathrm{tr}} = \boldsymbol{D}(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n^p), \quad \boldsymbol{I}_c = \boldsymbol{I} + u \boldsymbol{D} \boldsymbol{H}, \quad u = \frac{\Delta\gamma}{\varphi_{n+1}} \quad (13) $$

where $I$ is a fourth order unit tensor.

$I_c$ is the fourth order tensor that relates the final stress tensor $\boldsymbol{\sigma}_{n+1}$ to the trial stress tensor $\boldsymbol{\sigma}^{\mathrm{tr}}$

$$ \boldsymbol{\sigma}_{n+1} = \boldsymbol{I}_c^{-1} \boldsymbol{\sigma}^{\mathrm{tr}} \quad (14) $$

The Newton-Raphson iterative scheme is used for the local Equation (12). To preserve the quadratic rate of convergence that characterizes Newton's method, it is crucial to use the consistent tangent modulus. The derivation of the elastoplastic consistent tangent modulus is relatively lengthy but follows only usual application of consistent linearization concepts. The details of derivations are omitted and the final form of the elastoplastic consistent tangent modulus for the proposed constitutive model is as follows
$$ S = \begin{bmatrix} D = S^{-1} \\ \begin{bmatrix} \frac{1}{E_1} & -\frac{\nu_{12}}{E_1} & -\frac{\nu_{13}}{E_1} & 0 & 0 & 0 \\ \frac{1}{E_2} & -\frac{\nu_{23}}{E_2} & 0 & 0 & 0 & 0 \\ \frac{1}{E_3} & 0 & 0 & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{\frac{1}{G_{12}} & 0 & 0} \\ \multicolumn{6}{c}{\frac{1}{G_{13}} & 0} \\ \multicolumn{6}{c}{\frac{1}{G_{23}}} \end{bmatrix} \end{bmatrix}, \quad (3) $$

where $(E_1, E_2, E_3)$ are longitudinal and transverse Young's moduli, $(G_{12}, G_{13}, G_{23})$ are shear moduli and $(\nu_{12}, \nu_{13}, \nu_{23})$ are Poisson ratios. With the plane stress assumption, the elasticity matrix $D$ became

$$ D = \begin{bmatrix} H_{11} & H_{12} & 0 & 0 & 0 \\ H_{22} & 0 & 0 & 0 & 0 \\ 0 & G_{12} & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{G_{13}} \\ \multicolumn{6}{c}{G_{23}} \end{bmatrix} \quad (4) $$

where $G_{12} = G_{13}$. The components of $D$ can be written in terms of Engineering constants as follows:

$$ \left\{ \begin{array}{l} H_{11} = \frac{E_1}{a}, H_{22} = \frac{E_2}{a}, H_{12} = \frac{\nu_{12}E_2}{a} = \frac{\nu_{21}E_1}{a} \\ a = 1 - \nu_{12}\nu_{21} = 1 - \frac{\nu_{12}^2 E_2}{E_1} \end{array} \right. \quad (5) $$

Concerning the yield criterion, it is chosen as

$$ f = \varphi(\boldsymbol{\sigma}) - \sigma_p, \quad \sigma_p = \sigma_y + R \quad (6) $$

where $\sigma_y$ is the initial yield stress, $\varphi(\boldsymbol{\sigma})$ is the equivalent stress and $R$ is the drag stress in isotropic hardening. In order to express the anisotropic plasticity, the equivalent stress $\varphi(\boldsymbol{\sigma})$ can be described by using the Hill criterion [26]

$$ \varphi(\boldsymbol{\sigma}) = \|\boldsymbol{\sigma}\|_{H} = \sqrt{\boldsymbol{\sigma}^{T} \boldsymbol{H} \boldsymbol{\sigma}} \quad (7) $$

where $H$ is a fourth-order tensor given in a condensed form as [27-29].

$$ H = \begin{bmatrix} H+G & -H & -G & 0 & 0 & 0 \\ H+F & -F & 0 & 0 & 0 & 0 \\ F+G & 0 & 0 & 0 & 0 & 0 \\ \multicolumn{6}{c}{Sym} \\ \multicolumn{6}{c}{2N} \\ \multicolumn{6}{c}{2M} \\ \multicolumn{6}{c}{2L} \end{bmatrix} \quad (8) $$

$H, G, F, N, M$ and $L$ are Hill material constants to be identified. The evolution equations for the present problem are given as follows

$$ \dot{\boldsymbol{\epsilon}}^p = \dot{\gamma} \frac{\partial f}{\partial \boldsymbol{\sigma}} = \dot{\gamma} \boldsymbol{n}, \quad \boldsymbol{n} = \frac{1}{\varphi(\boldsymbol{\sigma})} \boldsymbol{H} \boldsymbol{\sigma} \quad (9) $$

$$ \dot{r} = -\dot{\gamma} \frac{\partial f}{\partial R} = \dot{\gamma} \quad (10) $$

where $\dot{\gamma}$ is the plastic multiplier which is consistent with the loading/unloading condition

$$ \dot{\gamma} \ge 0, \quad f \le 0, \quad \dot{\gamma} f = 0 \quad (11) $$

## 2.2 Numerical Integration

A fully implicit integration algorithm is used to solve the constitutive equations. The main advantage of the present local fully implicit integration algorithm is its strain driven characteristic and unconditional stability with a low number of equations to be solved: only one equation. The unknown of this equation is the plastic multiplier $\Delta\gamma$

$$ f(\Delta\gamma) = \varphi(\boldsymbol{\sigma}) - \sigma_p = 0, \quad \varphi = \left[ \boldsymbol{\sigma}^{\mathrm{tr}T} \boldsymbol{I}_c^{-T} \boldsymbol{H} \boldsymbol{I}_c^{-1} \boldsymbol{\sigma}^{\mathrm{tr}} \right]^{1/2} \quad (12) $$

$$ \boldsymbol{\sigma}^{\mathrm{tr}} = \boldsymbol{D}(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}_n^p), \quad \boldsymbol{I}_c = \boldsymbol{I} + u \boldsymbol{D} \boldsymbol{H}, \quad u = \frac{\Delta\gamma}{\varphi_{n+1}} \quad (13) $$

where $I$ is a fourth order unit tensor.

$I_c$ is the fourth order tensor that relates the final stress tensor $\boldsymbol{\sigma}_{n+1}$ to the trial stress tensor $\boldsymbol{\sigma}^{\mathrm{tr}}$

$$ \boldsymbol{\sigma}_{n+1} = \boldsymbol{I}_c^{-1} \boldsymbol{\sigma}^{\mathrm{tr}} \quad (14) $$

The Newton-Raphson iterative scheme is used for the local Equation (12). To preserve the quadratic rate of convergence that characterizes Newton's method, it is crucial to use the consistent tangent modulus. The derivation of the elastoplastic consistent tangent modulus is relatively lengthy but follows only usual application of consistent linearization concepts. The details of derivations are omitted and the final form of the elastoplastic consistent tangent modulus for the proposed constitutive model is as follows
$$ \left. \frac{d\sigma_{n+1}}{d\varepsilon_{n+1}} \right|_{n+1} = D^{ep} = D^* - \frac{1}{q} V_1 \otimes V_2 \quad (15) $$

$$ D^* = I_c^{-1} D, \quad V_1 = u_1 \varphi_{n+1} D^* n_{n+1}, \quad V_2 = D^{*T} n_{n+1} \quad (16) $$

$$ q = \sigma_p' + n_{n+1}^T V_1 \quad (17) $$

The numerical treatment of the proposed algorithm is implemented in ABAQUS using user interface material subroutines UMAT. More details on the formulation and numerical implementation of the present model can be found in Wali et al. [27] and Wali et al. [30].

## 2.3 Extensions to Viscoplasticity

As proposed by Simo et al. [31] and followed by Koubaa et al. [32], the loading/unloading condition, Equation (12) for rate independent plasticity, is rewritten, for rate dependent plasticity, in terms of the inverse function $g^{-1}(.)$ as

$$ f = g^{-1}\left(\frac{\Delta\gamma}{\Delta t}\right) \quad \text{for } f > 0 \Leftrightarrow \Delta\gamma > 0 \quad (18) $$

Combining Equation (18) with the consistency condition Equation (12) results in the modified Equation (19)

$$ h(\Delta\gamma) = \varphi(\sigma) - \left[ \sigma_p + g^{-1}\left(\frac{\Delta\gamma}{\Delta t}\right) \right] = 0 \quad (19) $$

For the linearization of Equation (12), and tangent operator, all to do is to replace

$$ (\sigma_p, \sigma_p') \text{ by } \left( \sigma_p + \bar{g} \left( \frac{\Delta\gamma}{\Delta t} \right), \sigma_p' + \frac{1}{\Delta t} \bar{g}' \left( \frac{\Delta\gamma}{\Delta t} \right) \right) \quad (20) $$

where $\bar{g}(.) = g^{-1}(.)$.

One can note that for practically all rate-dependent models through the choice of $g(.)$, the inverse $g^{-1}(.)$ is given in a closed form. In the present paper a modified Johnson-Cook (J.C.) law is used where the $g(.)$ function is given by

$$ g(x) = \dot{\varepsilon}_0 \exp \left[ \frac{x}{C \sigma_p} \right] \quad (21) $$

which leads to the following equation for the flow stress

$$ h = \varphi(\sigma) - \left[ \sigma_y + Q \left( 1 - \exp(-\beta \varepsilon_p) \right) + k \varepsilon_p \right] \left[ 1 - C \text{Log} \left( \frac{\dot{\varepsilon}}{\dot{\varepsilon}_0} \right) \right] = 0 \quad (22) $$

$k$ is the asymptotic hardening modulus, $Q$ and $\beta$ describe the non-linear part of the tensile curve when plasticity appears, and $C$ represents the strain rate

sensitivity coefficient. $\dot{\varepsilon}_0$ is the reference strain rate and $\dot{\varepsilon}$ represents the strain rate.

# 3 MODEL IDENTIFICATION AND VALIDATION IN QUASI-STATIC CONDITIONS

## 3.1 Anisotropic Elasticity

The tested material is a quasi-unidirectional flax fiber reinforced epoxy composite, based on so-called FUD 180 flax reinforcement supplied by LINEO NV Company. The quasi-unidirectional flax/epoxy composite is marked by an anisotropic behavior. Using tensile tests on samples with different fiber orientations (0°, 45° and 90°), the elastic properties of this composite are identified and given in Table 1.

Tensile tests at different strain rates performed by Poilane et al. [7] are explored in order to identify the behavior of flax fiber reinforced epoxy. The strain rates are $\dot{\varepsilon}_0 = 10^{-6} s^{-1}$; $\dot{\varepsilon} = 10^{-5} s^{-1}$; $\dot{\varepsilon} = 10^{-4} s^{-1}$; $\dot{\varepsilon} = 10^{-3} s^{-1}$ and the stress-strain curves are presented in Poilane et al. [7]. According to these tests, flax fiber reinforced composites are characterized by a viscoplastic behavior.

<table><caption>**Table 1** Elastic properties of FUD180 reinforced composite materials [8].</caption><thead><tr><th>E<sub>11</sub> (GPa)</th><th>E<sub>22</sub> = E<sub>33</sub> (GPa)</th><th>G<sub>12</sub> = G<sub>13</sub> (GPa)</th><th>G<sub>23</sub> (GPa)</th><th>ν<sub>12</sub> = ν<sub>13</sub></th><th>ν<sub>23</sub></th></tr></thead><tbody><tr><td>18.0±</td><td>4.7±</td><td>1.86±</td><td>1.66±</td><td>0.37±</td><td>0.40±</td></tr><tr><td>1.0</td><td>0.3</td><td>0.17</td><td>0.15</td><td>0.01</td><td>0.06</td></tr></tbody></table>

## 3.2 Hill's Material Parameters Identification

Tensile stress-strain curves obtained with composite samples at different fiber orientations (0°, 45°, and 90°) (Figure 1) are used to determine the Hill's material parameters, [26], Table 2. These material parameters are determined via equations that relate the yield stresses for different fiber orientations ($\sigma_{0°}$, $\sigma_{45°}$, $\sigma_{90°}$) to one another and with eventually the coefficient $r_0$ [26]. This coefficient $r_0$ is the ratio between the two transverse strains (strains among the width and the thickness). In this case, $r_0=1$.

<table><caption>**Table 2** Hill's material parameters.</caption><thead><tr><th>F</th><th>G</th><th>H</th><th>N</th><th>M</th><th>L</th></tr></thead><tbody><tr><td>1 / (<sup>σ<sub>90°</sub></sup>/<sub>σ<sub>0°</sub></sub>)<sup>2</sup> - H</td><td>1 / (1 + r<sub>0</sub>)</td><td>1 - G</td><td>1 / 2 ⋅ (<sup>4</sup>/<sub>(σ<sub>45°</sub>/σ<sub>0°</sub>)<sup>2</sup> - F - G</sub>)</td><td>N</td><td>N</td></tr><tr><td>2.51</td><td>0.50</td><td>0.50</td><td>1.97</td><td>1.97</td><td>1.97</td></tr></tbody></table>
$$ \left. \frac{d\sigma_{n+1}}{d\varepsilon_{n+1}} \right|_{n+1} = D^{ep} = D^* - \frac{1}{q} V_1 \otimes V_2 \quad (15) $$

$$ D^* = I_c^{-1} D, \quad V_1 = u_1 \varphi_{n+1} D^* n_{n+1}, \quad V_2 = D^{*T} n_{n+1} \quad (16) $$

$$ q = \sigma_p' + n_{n+1}^T V_1 \quad (17) $$

The numerical treatment of the proposed algorithm is implemented in ABAQUS using user interface material subroutines UMAT. More details on the formulation and numerical implementation of the present model can be found in Wali et al. [27] and Wali et al. [30].

## 2.3 Extensions to Viscoplasticity

As proposed by Simo et al. [31] and followed by Koubaa et al. [32], the loading/unloading condition, Equation (12) for rate independent plasticity, is rewritten, for rate dependent plasticity, in terms of the inverse function $g^{-1}(.)$ as

$$ f = g^{-1}\left(\frac{\Delta\gamma}{\Delta t}\right) \quad \text{for } f > 0 \Leftrightarrow \Delta\gamma > 0 \quad (18) $$

Combining Equation (18) with the consistency condition Equation (12) results in the modified Equation (19)

$$ h(\Delta\gamma) = \varphi(\sigma) - \left[ \sigma_p + g^{-1}\left(\frac{\Delta\gamma}{\Delta t}\right) \right] = 0 \quad (19) $$

For the linearization of Equation (12), and tangent operator, all to do is to replace

$$ (\sigma_p, \sigma_p') \text{ by } \left( \sigma_p + \bar{g}\left(\frac{\Delta\gamma}{\Delta t}\right), \sigma_p' + \frac{1}{\Delta t}\bar{g}'\left(\frac{\Delta\gamma}{\Delta t}\right) \right) \quad (20) $$

where $\bar{g}(.) = g^{-1}(.)$.

One can note that for practically all rate-dependent models through the choice of $g(.)$, the inverse $g^{-1}(.)$ is given in a closed form. In the present paper a modified Johnson-Cook (J.C.) law is used where the $g(.)$ function is given by

$$ g(x) = \dot{\varepsilon}_0 \exp \left[ \frac{x}{C \sigma_p} \right] \quad (21) $$

which leads to the following equation for the flow stress

$$ h = \varphi(\sigma) - \left[ \sigma_y + Q \left(1 - \exp(-\beta \varepsilon_p)\right) + k \varepsilon_p \right] \left[ 1 - C \text{Log} \left( \frac{\dot{\varepsilon}}{\dot{\varepsilon}_0} \right) \right] = 0 \quad (22) $$

$k$ is the asymptotic hardening modulus, $Q$ and $\beta$ describe the non-linear part of the tensile curve when plasticity appears, and $C$ represents the strain rate

sensitivity coefficient. $\dot{\varepsilon}_0$ is the reference strain rate and $\dot{\varepsilon}$ represents the strain rate.

# 3 MODEL IDENTIFICATION AND VALIDATION IN QUASI-STATIC CONDITIONS

## 3.1 Anisotropic Elasticity

The tested material is a quasi-unidirectional flax fiber reinforced epoxy composite, based on so-called FUD 180 flax reinforcement supplied by LINEO NV Company. The quasi-unidirectional flax/epoxy composite is marked by an anisotropic behavior. Using tensile tests on samples with different fiber orientations (0°, 45° and 90°), the elastic properties of this composite are identified and given in Table 1.

Tensile tests at different strain rates performed by Poilane et al. [7] are explored in order to identify the behavior of flax fiber reinforced epoxy. The strain rates are $\dot{\varepsilon}_0 = 10^{-6} s^{-1}$; $\dot{\varepsilon} = 10^{-5} s^{-1}$; $\dot{\varepsilon} = 10^{-4} s^{-1}$; $\dot{\varepsilon} = 10^{-3} s^{-1}$ and the stress-strain curves are presented in Poilane et al. [7]. According to these tests, flax fiber reinforced composites are characterized by a viscoplastic behavior.

<table><caption>**Table 1** Elastic properties of FUD180 reinforced composite materials [8].</caption><thead><tr><th>E<sub>11</sub> (GPa)</th><th>E<sub>22</sub> = E<sub>33</sub> (GPa)</th><th>G<sub>12</sub> = G<sub>13</sub> (GPa)</th><th>G<sub>23</sub> (GPa)</th><th>ν<sub>12</sub> = ν<sub>13</sub></th><th>ν<sub>23</sub></th></tr></thead><tbody><tr><td>18.0±</td><td>4.7±</td><td>1.86±</td><td>1.66±</td><td>0.37±</td><td>0.40±</td></tr><tr><td>1.0</td><td>0.3</td><td>0.17</td><td>0.15</td><td>0.01</td><td>0.06</td></tr></tbody></table>

## 3.2 Hill's Material Parameters Identification

Tensile stress-strain curves obtained with composite samples at different fiber orientations (0°, 45°, and 90°) (Figure 1) are used to determine the Hill's material parameters, [26], Table 2. These material parameters are determined via equations that relate the yield stresses for different fiber orientations ($\sigma_{0°}$, $\sigma_{45°}$, $\sigma_{90°}$) to one another and with eventually the coefficient $r_0$ [26]. This coefficient $r_0$ is the ratio between the two transverse strains (strains among the width and the thickness). In this case, $r_0=1$.

<table><caption>**Table 2** Hill's material parameters.</caption><thead><tr><th>F</th><th>G</th><th>H</th><th>N</th><th>M</th><th>L</th></tr></thead><tbody><tr><td>1 / (<sup>σ<sub>90°</sub></sup>/<sub>σ<sub>0°</sub></sub>)<sup>2</sup> - H</td><td>1 / (1 + r<sub>0</sub>)</td><td>1 - G</td><td>1 / 2 ⋅ (<sup>4</sup>/<sub>(σ<sub>45°</sub>/σ<sub>0°</sub>)<sup>2</sup></sub>) - F - G</td><td>N</td><td>N</td></tr><tr><td>2.51</td><td>0.50</td><td>0.50</td><td>1.97</td><td>1.97</td><td>1.97</td></tr></tbody></table>
<figure><img src="image_1.png" /><figcaption>**Figure 1** Stress-strain curves at different fiber orientations (0°, 45° and 90°) [8].</figcaption></figure>

### 3.3 Isotropic Hardening and Johnson Cook Parameters Identification

In this study, the drag stress *R* is expressed by an optimized Johnson Cook law defined in Equation (22). The identification of the modified J.C. law consists in determining the material parameters *Q*, *β*, *k*, and *C*.

This work is achieved in two steps. The first step consists in identifying the hardening parameters (*Q*, *β* and *k*) using the least square method for the stress-strain curve at the reference strain rate (ε̇₀ = 10⁻⁶ s⁻¹) (Figure 2). It is observed that the initiation of plastic deformation is not well modeled. However, J.C. model becomes completely effective when the flow is stabilized.

In the second step, the coefficient *C* is determined by inserting the values of the hardening parameters and those of the strain rates in Equation (22). The coefficient *C* is the slope of the linear trend curve obtained by plotting (σ<sub>exp</sub>/σ<sub>0</sub> - 1) as a function of the log of the deformation velocity (Log (ε̇/ε̇₀)). The identified material parameters of the optimized Johnson Cook law are given in Table 3. Finally, the model is validated with the experimental curves at different strain rates as presented in Figure 3.

<table><caption>**Table 3** Material parameters of the optimized Johnson Cook law.</caption><thead><tr><td>σ<sub>y</sub> (MPa)</td><td>Q (MPa)</td><td>β</td><td>k (MPa)</td><td>C</td></tr></thead><tbody><tr><td>33</td><td>6.015</td><td>1210</td><td>20684</td><td>0.152</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 2** Validation of the stress parameters at the reference strain rate (dashed line) with experimental curve at ε̇₀ = 10⁻⁶ s⁻¹ (orientation 0°).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 3** Validation of the model (dashed lines) with the experimental curves [7] (orientation 0°).</figcaption></figure>

### 3.4 Unidirectional Tensile Test Verification

Tensile test simulations using the J.C. model implemented in a finite element code Abaqus/Umat is accomplished in order to validate the model with experimental tensile test results. Finite element tensile tests are simulated under plane stress conditions. Figure 4 presents the initial configuration including the specimen shape and dimensions used for the numerical tensile simulations as well as the deformed configurations for different fiber orientations (θ = 0°, 45° and 90°) after tensile tests.

For a fiber orientation of 45° (Figure 4c), the deformed shape of tensile specimen is inclined. This phenomenon can probably be due to the anisotropy of the flax/epoxy composites. Tensile loading of the composite with fibers oriented at 45° may cause fiber/matrix shear which can
<figure><img src="image_3.png" /><figcaption>**Figure 1** Stress-strain curves at different fiber orientations (0°, 45° and 90°) [8].</figcaption></figure>

### 3.3 Isotropic Hardening and Johnson Cook Parameters Identification

In this study, the drag stress *R* is expressed by an optimized Johnson Cook law defined in Equation (22). The identification of the modified J.C. law consists in determining the material parameters *Q*, *β*, *k*, and *C*.

This work is achieved in two steps. The first step consists in identifying the hardening parameters (*Q*, *β* and *k*) using the least square method for the stress-strain curve at the reference strain rate (ε̇₀ = 10⁻⁶ s⁻¹) (Figure 2). It is observed that the initiation of plastic deformation is not well modeled. However, J.C. model becomes completely effective when the flow is stabilized.

In the second step, the coefficient *C* is determined by inserting the values of the hardening parameters and those of the strain rates in Equation (22). The coefficient *C* is the slope of the linear trend curve obtained by plotting (σ<sub>exp</sub>/σ<sub>0</sub> - 1) as a function of the log of the deformation velocity (Log (ε̇/ε̇₀)). The identified material parameters of the optimized Johnson Cook law are given in Table 3. Finally, the model is validated with the experimental curves at different strain rates as presented in Figure 3.

<table><caption>**Table 3** Material parameters of the optimized Johnson Cook law.</caption><thead><tr><td>σ<sub>y</sub> (MPa)</td><td>Q (MPa)</td><td>β</td><td>k (MPa)</td><td>C</td></tr></thead><tbody><tr><td>33</td><td>6.015</td><td>1210</td><td>20684</td><td>0.152</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2** Validation of the stress parameters at the reference strain rate (dashed line) with experimental curve at ε̇₀ = 10⁻⁶ s⁻¹ (orientation 0°).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3** Validation of the model (dashed lines) with the experimental curves [7] (orientation 0°).</figcaption></figure>

### 3.4 Unidirectional Tensile Test Verification

Tensile test simulations using the J.C. model implemented in a finite element code Abaqus/Umat is accomplished in order to validate the model with experimental tensile test results. Finite element tensile tests are simulated under plane stress conditions. Figure 4 presents the initial configuration including the specimen shape and dimensions used for the numerical tensile simulations as well as the deformed configurations for different fiber orientations (θ = 0°, 45° and 90°) after tensile tests.

For a fiber orientation of 45° (Figure 4c), the deformed shape of tensile specimen is inclined. This phenomenon can probably be due to the anisotropy of the flax/epoxy composites. Tensile loading of the composite with fibers oriented at 45° may cause fiber/matrix shear which can
also explain the inclination observed in the specimen. Figure 5 depicts a comparison between finite element and experimental tensile results. It is shown that for fiber orientation of 0°, the finite element tensile test simulations are in good agreement with the experimental results. As for fiber orientation of 45° and 90°, the numerical tensile curves are valid in the elastic domain and at the beginning of the plasticity. This can be explained by the fact that the present formulation does not take into account the damage phenomena. One other explanation could be that for fiber orientation of 45° and 90°, the viscous effect is more pronounced as the matrix play a crucial role in those directions. Including damage in the formulation as in [33-35] would be an interesting topic in predicting the exact behavior of the quasi-unidirectional flax (FUD180)/epoxy composites. It would be also useful to take into account of the viscosity at several characteristic scales [36, 37] or the coupling between damage and viscoelasticity as in [38].

<figure><img src="image_4.png" /><figcaption>**Figure 4** Unidirectional tensile tests: Finite element simulation for different fiber orientation ((a): Initial configuration; (b): $\theta = 0°$; (c): $\theta = 45°$; (d): $\theta = 90°$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5** Uniaxial tensile tests at 0°, 45° and 90°, present finite element results versus experimental results from Cherif et al. [8].</figcaption></figure>

## 3.5 Three-Point Bending Test

The purpose of this work is to validate the identified constitutive model, developed in Sub-sections 3.1, 3.2 and 3.3 with tests that have not served for the identification. The model is implemented in an Abaqus/Umat in order to reproduce acceptable results in the case of three-point bending test. An initially specimen shell supported on ends is subjected to a transverse displacement in the middle. The specimen length is L=80 mm, the width is w=10 mm and the thickness is $e = 2 mm$ (Figure 6). The material elastic properties are given in Table 1. Hill's material parameters are given in Table 2. As for the material parameters of J.C. law are brought up in Table 3.

The material investigated in the experimental and numerical study is flax fiber reinforced epoxy composite. The same material is used for tensile tests that have served for the identification of the J.C. model. The reinforcement is a quasi-unidirectional flax ply FUD 180 (weight 180 g/m²) manufactured by LINEO NV Company. The matrix is an epoxy system which consists of a combination of a resin (LY 1546) and a hardener (Aradur 3787) supplied by HUNTSMAN Corporation. The flax reinforcements are dried for one hour at 110°C in order to eliminate water molecules. The composite is manufactured using compression molding. The reinforcement is impregnated manually then kept under compression at ambient temperature for two hours before curing at 100°C for two hours. The samples are cut using Laser machine Trotec speedy 400.

Three-point bending tests have been carried out using an MTS criterion tensile testing machine Model 43. The distance between the two cylinders (diameter 25 mm) for fixing the test piece is 70 mm. Three-point bending tests are carried out in air at 20°C.
also explain the inclination observed in the specimen. Figure 5 depicts a comparison between finite element and experimental tensile results. It is shown that for fiber orientation of 0°, the finite element tensile test simulations are in good agreement with the experimental results. As for fiber orientation of 45° and 90°, the numerical tensile curves are valid in the elastic domain and at the beginning of the plasticity. This can be explained by the fact that the present formulation does not take into account the damage phenomena. One other explanation could be that for fiber orientation of 45° and 90°, the viscous effect is more pronounced as the matrix play a crucial role in those directions. Including damage in the formulation as in [33-35] would be an interesting topic in predicting the exact behavior of the quasi-unidirectional flax (FUD180)/epoxy composites. It would be also useful to take into account of the viscosity at several characteristic scales [36, 37] or the coupling between damage and viscoelasticity as in [38].

<figure><img src="image_6.png" /><figcaption>**Figure 4** Unidirectional tensile tests: Finite element simulation for different fiber orientation ((a): Initial configuration; (b): $\theta = 0°$; (c): $\theta = 45°$; (d): $\theta = 90°$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5** Uniaxial tensile tests at 0°, 45° and 90°, present finite element results versus experimental results from Cherif et al. [8].</figcaption></figure>

## 3.5 Three-Point Bending Test

The purpose of this work is to validate the identified constitutive model, developed in Sub-sections 3.1, 3.2 and 3.3 with tests that have not served for the identification. The model is implemented in an Abaqus/Umat in order to reproduce acceptable results in the case of three-point bending test. An initially specimen shell supported on ends is subjected to a transverse displacement in the middle. The specimen length is L=80 mm, the width is w=10 mm and the thickness is e = 2 mm (Figure 6). The material elastic properties are given in Table 1. Hill's material parameters are given in Table 2. As for the material parameters of J.C. law are brought up in Table 3.

The material investigated in the experimental and numerical study is flax fiber reinforced epoxy composite. The same material is used for tensile tests that have served for the identification of the J.C. model. The reinforcement is a quasi-unidirectional flax ply FUD 180 (weight 180 g/m²) manufactured by LINEO NV Company. The matrix is an epoxy system which consists of a combination of a resin (LY 1546) and a hardener (Aradur 3787) supplied by HUNTSMAN Corporation. The flax reinforcements are dried for one hour at 110°C in order to eliminate water molecules. The composite is manufactured using compression molding. The reinforcement is impregnated manually then kept under compression at ambient temperature for two hours before curing at 100°C for two hours. The samples are cut using Laser machine Trotec speedy 400.

Three-point bending tests have been carried out using an MTS criterion tensile testing machine Model 43. The distance between the two cylinders (diameter 25 mm) for fixing the test piece is 70 mm. Three-point bending tests are carried out in air at 20°C.
The numerical model of the three-point bending test is shown in Figure 6. The numerical study is conducted with the commercial software ABAQUS. The specimen is modeled with the standard structural shell elements S4 using 2400 elements and 10 integration points in thickness. It is a quadrilateral 4-nodes element with three rotational and three translational degrees of freedom per node. The contacts punch/specimen and specimen/dies are referred to as “hard contact”. The friction coefficient is 0.2. A convergence study indicates that the use 5×40 nodes are appropriate. This problem is solved numerically using the UMAT implemented in ABAQUS/implicit.

Figure 7 presents the initial and the contour plots of the transverse displacement of the final real deformed configuration for the flax fiber reinforced epoxy specimen. A displacement of 10 mm is imposed. The maximum displacement is observed in the middle of the tested specimen which is in accordance with the experimental results.

A comparison of the induced load evolution in the middle of the specimen versus the imposed displacement between the experimental and finite element simulation for a strain rate in the range of 10⁻⁴ s⁻¹ at different fiber orientations is shown in Figures 8, 9 and 10.

<figure><img src="image_6.png" /><figcaption>**Figure 6** Numerical model of three points bending test.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7** Three points bending test. Z-displacement on the final real deformed configuration θ = 0°.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8** Load versus displacement responses for three points bending test for fiber orientation 0°.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 9** Load versus displacement responses for three points bending test for fiber orientation 45°.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10** Load versus displacement responses for three points bending test for fiber orientation 90°.</figcaption></figure>

Numerical and experimental results for fiber orientation 0°, 45° and 90° are in good agreement especially in the first stage, where the proposed model is applicable. These results illustrate the good performance of the proposed model to predict the mechanical behavior of the studied material. Conversely to three-point bending
The numerical model of the three-point bending test is shown in Figure 6. The numerical study is conducted with the commercial software ABAQUS. The specimen is modeled with the standard structural shell elements S4 using 2400 elements and 10 integration points in thickness. It is a quadrilateral 4-nodes element with three rotational and three translational degrees of freedom per node. The contacts punch/specimen and specimen/dies are referred to as “hard contact”. The friction coefficient is 0.2. A convergence study indicates that the use 5×40 nodes are appropriate. This problem is solved numerically using the UMAT implemented in ABAQUS/implicit.

Figure 7 presents the initial and the contour plots of the transverse displacement of the final real deformed configuration for the flax fiber reinforced epoxy specimen. A displacement of 10 mm is imposed. The maximum displacement is observed in the middle of the tested specimen which is in accordance with the experimental results.

A comparison of the induced load evolution in the middle of the specimen versus the imposed displacement between the experimental and finite element simulation for a strain rate in the range of 10⁻⁴ s⁻¹ at different fiber orientations is shown in Figures 8, 9 and 10.

<figure><img src="image_8.png" /><figcaption>**Figure 6** Numerical model of three points bending test.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7** Three points bending test. Z-displacement on the final real deformed configuration θ = 0°.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8** Load versus displacement responses for three points bending test for fiber orientation 0°.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9** Load versus displacement responses for three points bending test for fiber orientation 45°.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10** Load versus displacement responses for three points bending test for fiber orientation 90°.</figcaption></figure>

Numerical and experimental results for fiber orientation 0°, 45° and 90° are in good agreement especially in the first stage, where the proposed model is applicable. These results illustrate the good performance of the proposed model to predict the mechanical behavior of the studied material. Conversely to three-point bending
behavior, the tensile behavior at 0° is better reproduced than the other orientations. Experimentally, it is under bending at 0° that delamination is observed the most. In the second stage, there is a discrepancy between the numerical and experimental results. As the proposed model does not take into consideration damage phenomena, it can be an eventual cause of this divergence. Local evolution at the punch and dies levels can cause this divergence as well.

# 4 MODELING OF LOW VELOCITY IMPACT TESTS

## 4.1 Finite Element Model Description

In this section, falling weight impact tests on a flax/epoxy composite circular plate are simulated in order to test the model and to deal with impact phenomena encountered in the automotive and aeronautic fields. Impact is imparted perpendicularly to the plane of the laminated plate, using a falling weight of mass 54.4 g with a hemispherical nosed cylinder of 6.35 mm diameter. The purpose of using small mass with a relatively high velocity of impact (from 10 m/s to 20 m/s) is to reproduce the conditions of severe strikes during service (an example could be the dropping of a small rock on an exterior automotive panel). The projectile hits the flax/epoxy circular plate of 60 mm radius with an initial velocity $V_0$ [32, 39, 40]. Fibers are oriented in the X direction. The circular plate is clamped at its borders and the projectile is subjected only to a movement in the Z direction (Figure 11). A rigid body is assumed for the impactor with property assigned. The contact between the projectile and the plate is referred to as “hard” contact. The friction coefficient between the projectile and the plate is 0.2. The plate is modeled with the standard structural shell elements S4. The plate is meshed using 38350 S4 elements and 5 integration points in thickness. In order to obtain more accurate results, the meshing is refined in the contact zone of the circular plate and of the impactor. The mechanical properties of the plate are illustrated in Tables 1, 2 and 3.

Figure 12 illustrates the influence of the C parameter in the Johnson-Cook model on the impact force evolution for $m_p=54.4$ g, $R_p=6.35$ mm, $e=1$ mm and $V_0=15$ m/s. The addition of the viscosity component in the (J.C.) law decreases the impact force which can be explained by the fact that a part of the impact energy is absorbed by the viscosity [40, 41].

<figure><img src="image_11.png" /><figcaption>Figure 11 Hemi-spherical projectile impacted circular plate.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12 Force versus time responses for $m_p=54.4$ g, $R_p=6.35$ mm, $e=1$ mm and $V_0=15$ m/s.</figcaption></figure>

## 4.2 Parametric Study

A parametric study is carried out in order to determine the influence of the plate thickness, the viscous component in the model, and the projectile's mass, velocity and radius on the response of the plate. The projectile hits the plate with an initial velocity $V_0$. Impact test can be divided into two phases. The first one consists in a loading phase where the projectile contacts the plate until its velocity and thus its energy reaches zero. This phase can be divided in two phases namely the elastic loading and the plastic loading. Then, due to the absorbed energy, the plate rebounds and pushes back the impactor. This second phase is referred to as unloading or recovery phase [29, 41].

* Influence of the plate thickness

As a first step, a study of the plate thickness effect on the plate reaction is presented in Figures 13, 14 and 15. The evolution of the plate center displacement with impact time is presented in Figure 13 for plate
behavior, the tensile behavior at 0° is better reproduced than the other orientations. Experimentally, it is under bending at 0° that delamination is observed the most. In the second stage, there is a discrepancy between the numerical and experimental results. As the proposed model does not take into consideration damage phenomena, it can be an eventual cause of this divergence. Local evolution at the punch and dies levels can cause this divergence as well.

# 4 MODELING OF LOW VELOCITY IMPACT TESTS

## 4.1 Finite Element Model Description

In this section, falling weight impact tests on a flax/epoxy composite circular plate are simulated in order to test the model and to deal with impact phenomena encountered in the automotive and aeronautic fields. Impact is imparted perpendicularly to the plane of the laminated plate, using a falling weight of mass 54.4 g with a hemispherical nosed cylinder of 6.35 mm diameter. The purpose of using small mass with a relatively high velocity of impact (from 10 m/s to 20 m/s) is to reproduce the conditions of severe strikes during service (an example could be the dropping of a small rock on an exterior automotive panel). The projectile hits the flax/epoxy circular plate of 60 mm radius with an initial velocity $V_0$ [32, 39, 40]. Fibers are oriented in the X direction. The circular plate is clamped at its borders and the projectile is subjected only to a movement in the Z direction (Figure 11). A rigid body is assumed for the impactor with property assigned. The contact between the projectile and the plate is referred to as “hard” contact. The friction coefficient between the projectile and the plate is 0.2. The plate is modeled with the standard structural shell elements S4. The plate is meshed using 38350 S4 elements and 5 integration points in thickness. In order to obtain more accurate results, the meshing is refined in the contact zone of the circular plate and of the impactor. The mechanical properties of the plate are illustrated in Tables 1, 2 and 3.

Figure 12 illustrates the influence of the C parameter in the Johnson-Cook model on the impact force evolution for $m_p=54.4$ g, $R_p=6.35$ mm, $e=1$ mm and $V_0=15$ m/s. The addition of the viscosity component in the (J.C.) law decreases the impact force which can be explained by the fact that a part of the impact energy is absorbed by the viscosity [40, 41].

<figure><img src="image_13.png" /><figcaption>Figure 11 Hemi-spherical projectile impacted circular plate.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12 Force versus time responses for $m_p=54.4$ g, $R_p=6.35$ mm, $e=1$ mm and $V_0=15$ m/s.</figcaption></figure>

## 4.2 Parametric Study

A parametric study is carried out in order to determine the influence of the plate thickness, the viscous component in the model, and the projectile's mass, velocity and radius on the response of the plate. The projectile hits the plate with an initial velocity $V_0$. Impact test can be divided into two phases. The first one consists in a loading phase where the projectile contacts the plate until its velocity and thus its energy reaches zero. This phase can be divided in two phases namely the elastic loading and the plastic loading. Then, due to the absorbed energy, the plate rebounds and pushes back the impactor. This second phase is referred to as unloading or recovery phase [29, 41].

### * Influence of the plate thickness

As a first step, a study of the plate thickness effect on the plate reaction is presented in Figures 13, 14 and 15. The evolution of the plate center displacement with impact time is presented in Figure 13 for plate
thicknesses equal to e=1 mm, e=2 mm and e=3 mm. According to the displacement response shown in Figure 13, it can be noted that the central deflection keeps on decreasing when the thickness increases. The impact force evolution in time predicted by the present numerical model for e=1 mm, e=2 mm and e=3 mm are plotted in Figure 14. The increase of the plate thickness causes an increase in the values of the maximal impact forces and a great decrease of the impact time. In fact, when the plate is thicker, the bending stiffness is higher which may explain the force increase and the deflection decrease. A reduction of the plate thickness increases the plate deformation (cf. Figure 21). Therefore, the specimen would need less time to find its equilibrium position in the unloading phase. This may be a possible cause of impact time decrease when the plate thickness increases. These results are in accordance with those shown in Figure 15 which illustrates that impact velocity reaches zero in less time as plate thickness increases. The contact force and impact velocity obtained from both numerical analyses indicate that the impact force will reach its peak when the velocity decreases to zero.

* **Influence of the projectile mass**

The analysis of the impact force evolution through time in Figure 16 indicates that the impactor mass increase augments the peak force and the impact time. As the projectile gets heavier, the impact energy is intensified which can explain the increase of the impact force. It may also justify the impact time increase since it becomes harder for the plate to push back the projectile. Figure 17 depicts the projectile mass effect on the evolution of the impact velocity in time. Velocity evolution is induced by the rigidity of the plate/impactor system and this one does not seem to change too much when the mass changes. The variation of the impactor's mass ($m_p$=23.5 g, $m_p$=35.5 g and $m_p$=54.4 g) are in accordance with the results obtained in the study of the plate thickness effect. Indeed, the impact velocity reaches zero at the impact force peak.

* **Influence of the projectile velocity**

Figure 18 depicts the effect of the projectile velocity on the impact force evolution through impact time. It is shown that the increase of velocity impact leads to an intensification of the impact energy and thus an increase in the impact force plate response.

Figure 19 illustrates the relationship between the maximum impact force and the initial impactor velocity for diverse mass of impactor. It is noticed that the maximum impact force increases with the increase of the initial impactor velocity, and can be expressed approximately in power relationship as follows,

$$F_{\max} = C_p V_0^{n_p} \qquad (23)$$

where, $C_p$ is a coefficient related to the impactor mass. The coefficients $C_p$ are given in Table 4. It can be noted that the projectile mass increase induces a drop in the values of the coefficients $C_p$ and $n_p$.

<table><caption>**Table 4** $C_p$ coefficient related to the projectile mass</caption><tr><td>$m_p$ (g)</td><td>23.5</td><td>35.55</td><td>54.4</td></tr><tr><td>$C_p$</td><td>0.0424</td><td>0.0365</td><td>0.0259</td></tr><tr><td>$n_p$</td><td>1.447</td><td>1.38</td><td>1.374</td></tr></table>

* **Anisotropy illustration**

Figure 20 shows deformation profiles in longitudinal and transverse directions of clamped circular plate for e=1 mm, where the fibers are oriented in the X direction. The difference observed between the X and the Y directions can be probably due to the orthotropic behavior of the flax/epoxy composite. These results are confirmed by the deformed configurations of flax/epoxy composite portrayed in Figure 21 at the time 0.88 ms under the initial velocity $V_0$=15 m/s for e=1 mm, e=2 mm and e=3 mm.

The anisotropy presented in composite materials affect greatly the low velocity impact behavior and is accompanied with a great variation of the maximum impact force and the center plate displacement of the impacted circular plate.

The contact impact model proposed in this paper can predict certain characteristics of the contact impact process. It is also able to predict the history of the impact force in the case of a projectile striking a circular plate.

<figure><img src="image_13.png" /><figcaption>**Figure 13** Displacement of the plate center versus time responses for $m_p$=54.4 g and $R_p$=6.35 mm and $V_0$=15 m/s.</figcaption></figure>
thicknesses equal to e=1 mm, e=2 mm and e=3 mm. According to the displacement response shown in Figure 13, it can be noted that the central deflection keeps on decreasing when the thickness increases. The impact force evolution in time predicted by the present numerical model for e=1 mm, e=2 mm and e=3 mm are plotted in Figure 14. The increase of the plate thickness causes an increase in the values of the maximal impact forces and a great decrease of the impact time. In fact, when the plate is thicker, the bending stiffness is higher which may explain the force increase and the deflection decrease. A reduction of the plate thickness increases the plate deformation (cf. Figure 21). Therefore, the specimen would need less time to find its equilibrium position in the unloading phase. This may be a possible cause of impact time decrease when the plate thickness increases. These results are in accordance with those shown in Figure 15 which illustrates that impact velocity reaches zero in less time as plate thickness increases. The contact force and impact velocity obtained from both numerical analyses indicate that the impact force will reach its peak when the velocity decreases to zero.

* **Influence of the projectile mass**

The analysis of the impact force evolution through time in Figure 16 indicates that the impactor mass increase augments the peak force and the impact time. As the projectile gets heavier, the impact energy is intensified which can explain the increase of the impact force. It may also justify the impact time increase since it becomes harder for the plate to push back the projectile. Figure 17 depicts the projectile mass effect on the evolution of the impact velocity in time. Velocity evolution is induced by the rigidity of the plate/impactor system and this one does not seem to change too much when the mass changes. The variation of the impactor's mass ($m_p$=23.5 g, $m_p$=35.5 g and $m_p$=54.4 g) are in accordance with the results obtained in the study of the plate thickness effect. Indeed, the impact velocity reaches zero at the impact force peak.

* **Influence of the projectile velocity**

Figure 18 depicts the effect of the projectile velocity on the impact force evolution through impact time. It is shown that the increase of velocity impact leads to an intensification of the impact energy and thus an increase in the impact force plate response.

Figure 19 illustrates the relationship between the maximum impact force and the initial impactor velocity for diverse mass of impactor. It is noticed that the maximum impact force increases with the increase of the initial impactor velocity, and can be expressed approximately in power relationship as follows,

$$F_{\max} = C_p V_0^{n_p} \qquad (23)$$

where, $C_p$ is a coefficient related to the impactor mass. The coefficients $C_p$ are given in Table 4. It can be noted that the projectile mass increase induces a drop in the values of the coefficients $C_p$ and $n_p$.

<table><caption>**Table 4** $C_p$ coefficient related to the projectile mass</caption><tr><td>$m_p$ (g)</td><td>23.5</td><td>35.55</td><td>54.4</td></tr><tr><td>$C_p$</td><td>0.0424</td><td>0.0365</td><td>0.0259</td></tr><tr><td>$n_p$</td><td>1.447</td><td>1.38</td><td>1.374</td></tr></table>

* **Anisotropy illustration**

Figure 20 shows deformation profiles in longitudinal and transverse directions of clamped circular plate for e=1 mm, where the fibers are oriented in the X direction. The difference observed between the X and the Y directions can be probably due to the orthotropic behavior of the flax/epoxy composite. These results are confirmed by the deformed configurations of flax/epoxy composite portrayed in Figure 21 at the time 0.88 ms under the initial velocity $V_0$=15 m/s for e=1 mm, e=2 mm and e=3 mm.

The anisotropy presented in composite materials affect greatly the low velocity impact behavior and is accompanied with a great variation of the maximum impact force and the center plate displacement of the impacted circular plate.

The contact impact model proposed in this paper can predict certain characteristics of the contact impact process. It is also able to predict the history of the impact force in the case of a projectile striking a circular plate.

<figure><img src="image_15.png" /><figcaption>**Figure 13** Displacement of the plate center versus time responses for $m_p$=54.4 g and $R_p$=6.35 mm and $V_0$=15 m/s.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 14** Impact force versus time responses for $m_p=54.4$ g and $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 15** Impact velocity versus time responses for $m_p=54.4$ g and $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 16** Impact force versus time responses for $e=1$ mm, $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 17** Impact velocity versus time responses for $e=1$ mm, $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 18** Impact force versus time responses for $e=1$ mm, $R_p=6.35$ mm and $m_p=54.4$ m/s.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 19** Maximum impact force versus initial impactor velocity.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14** Impact force versus time responses for $m_p=54.4$ g and $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15** Impact velocity versus time responses for $m_p=54.4$ g and $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16** Impact force versus time responses for $e=1$ mm, $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17** Impact velocity versus time responses for $e=1$ mm, $R_p=6.35$ mm and $V_0=15$ m/s.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18** Impact force versus time responses for $e=1$ mm, $R_p=6.35$ mm and $m_p=54.4$ m/s.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19** Maximum impact force versus initial impactor velocity.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 20** Deformation profiles at 0.88 ms impact time along the rolling and transverse directions for *m*<sub>p</sub>=54.4 g, *V*<sub>0</sub>=15 m/s and *e*=1 mm.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 21** Hill's equivalent strain at *t*=0.88 ms for *m*<sub>p</sub>=54.4 g and *V*<sub>0</sub>=15 m/s.</figcaption></figure>

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

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

**5 CONCLUSION**

Experimental tests have shown that flax fiber reinforced
polymer composites are characterized by nonlinear
visco-elastoplastic behavior. In this paper, a three
dimensional formulation with orthotropic elasticity and
orthotropic plasticity using Hill criterion is accomplished
to model the behavior of a quasi-unidirectional flax fiber
reinforced epoxy. The isotropic hardening and Johnson
Cook parameters are identified from unidirectional
tensile tests at different strain rates. A numerical
integration of the constitutive equations is implemented
in a user-defined material subroutines ABAQUS/UMAT.
The model validation is achieved via comparison
between experimental and finite element tensile tests
simulations for different fiber orientations. To validate
the numerical implementation, experimental and
numerical three-point bending tests are carried out. A
good agreement is found between the experimental and
finite element simulation results for different fiber
orientations. The discrepancy observed at high strain is
explained by the fact that damage is not taken into
account in this present model. Finally, a numerical
parametric study on low velocity impact of a flax/epoxy
composite circular plate is studied. A parametric study,
such as the plate thickness, the projectile mass, velocity
and radius are investigated.

In future work, in order to reproduce the complex
behavior of the flax/epoxy composites, enhancing the
proposed model and thus implementation are
considered by taking into account progressive damage
phenomena.
<figure><img src="image_22.png" /><figcaption>**Figure 20** Deformation profiles at 0.88 ms impact time along the rolling and transverse directions for *m*<sub>p</sub>=54.4 g, *V*<sub>0</sub>=15 m/s and *e*=1 mm.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21** Hill's equivalent strain at *t*=0.88 ms for *m*<sub>p</sub>=54.4 g and *V*<sub>0</sub>=15 m/s.</figcaption></figure>

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

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

## 5 CONCLUSION

Experimental tests have shown that flax fiber reinforced polymer composites are characterized by nonlinear visco-elastoplastic behavior. In this paper, a three dimensional formulation with orthotropic elasticity and orthotropic plasticity using Hill criterion is accomplished to model the behavior of a quasi-unidirectional flax fiber reinforced epoxy. The isotropic hardening and Johnson Cook parameters are identified from unidirectional tensile tests at different strain rates. A numerical integration of the constitutive equations is implemented in a user-defined material subroutines ABAQUS/UMAT. The model validation is achieved via comparison between experimental and finite element tensile tests simulations for different fiber orientations. To validate the numerical implementation, experimental and numerical three-point bending tests are carried out. A good agreement is found between the experimental and finite element simulation results for different fiber orientations. The discrepancy observed at high strain is explained by the fact that damage is not taken into account in this present model. Finally, a numerical parametric study on low velocity impact of a flax/epoxy composite circular plate is studied. A parametric study, such as the plate thickness, the projectile mass, velocity and radius are investigated.

In future work, in order to reproduce the complex behavior of the flax/epoxy composites, enhancing the proposed model and thus implementation are considered by taking into account progressive damage phenomena.
REFERENCES

1. H.L. Bos, The potential of flax fibres as reinforcement for composite materials (Ph.D. Thesis). Eindhoven: Technische Universiteit Eindhoven (2004).

2. C. Bailey, Analysis of the flax fibres tensile behavior and analysis of the tensile stiffness increase. *Composite Part A* **33**, 939-948 (2002).

3. C. Gourier, A. Le Duigou, A. Bourmaud and C. Bailey, Mechanical analysis of elementary flax fibre tensile properties after different thermal cycles. *Composite Part A* **64**, 159-166 (2014).

4. M. Fehri, R.R. Ragueh, A. Vivet, F. Dammak and M. Haddar, Improvement of natural fiber composite materials by carbon fibers. *Journal of Renewable Materials* **5**(1), 38-47 (2017).

5. J. Andersons, J. Modniks and E. Sparnins, Modeling the nonlinear deformation of flax-fiber-reinforced polymer matrix laminates in active loading. *Journal of Reinforced Plastics and Composites* **34**(3), 248-256 (2015).

6. A. Rubio-López, T. Hoang and C. Santiuste, Constitutive model to predict the viscoplastic behaviour of natural fibres based composites. *Journal of Composite Structures* **155**, 8-18 (2016).

7. C. Poilane, Z.E. Cherif, F. Richard, A. Vivet, B. Ben Doudou and J. Chen, Polymer reinforced by flax fibres as a viscoelastoplastic material. *Journal of Composite Structures* **112**, 100-112 (2014).

8. Z.E. Cherif, C. Poilâne, A. Vivet, B. Ben Doudou and J. Chen. About optimal architecture of plant fibre textile composite for mechanical and sorption properties. *Journal of Composite Structures* **140**, 240-251 (2016).

9. M.J. John and S. Thomas, Biofibres and biocomposites. *Carbohydr Polym* **71**, 343-364 (2008).

10. A.K. Bledzki and A. Jaszkiewicz, Mechanical performance of biocomposites based on PLA and PHBV reinforced with natural fibres-a comparative study to PP. *Composites Science and Technology* **70**, 1687-1696 (2010).

11. H. Ku, Y.M. Cheng, C. Snook and D. Baddeley, Drop weight impact test fracture of vinyl ester composites: micrographs of pilot study. *Journal of Composite Materials* **39**, 1607-1620 (2005).

12. C. Santiuste, S. Sanchez-Saez and E. Barbero, Residual flexural strength after low velocity impact in glass/polyester composite beams. *Journal of Composite Structures* **92**, 25-30 (2010).

13. C. Santiuste, S. Sanchez-Saez and E. Barbero, A comparison of progressive-failure criteria in the prediction of the dynamic bending failure of composite laminated beams. *Journal of Composite Structures* **92**, 2406-2414 (2010).

14. I. Ivañez, C. Santiuste and S. Sanchez-Saez, FEM analysis of dynamic flexural behavior of composite sandwich beams with foam core. *Journal of Composite Structures* **92**, 2285-2291 (2010).

15. A.R. Setoodeh, P. Malekzadeh and K. Nikbin, Low velocity impact analysis of laminated composite plates using a 3D elasticity based layerwise FEM. *Materials & Design* **30**, 3795-3801 (2009).

16. S.H. Song, Y.S. Byun, T.W. Ku, W.J. Song, J. Kim and B.S. Kang, Experimental and numerical investigation on impact performance of carbon reinforced aluminum laminates. *Journal of Materials Science & Technology* **26**, 327-332 (2010).

17. S.M.R. Khalili, M. Soroush, A. Davar and O. Rahmani, Finite element modeling of low-velocity impact on laminated composite plates and cylindrical shells. *Journal of Composite Structures* **93**, 1363-1375 (2011).

18. D. Zhang, Y. Sun, L. Chen and N. Pan, A comparative study on low-velocity impact response of fabric composite laminates. *Materials & Design* **50**, 750-756 (2013).

19. L. Ferrante, F. Sarasini, J. Tirillò, L. Lampani, T. Valente and P. Gaudenzi, Low velocity impact response of basalt-aluminium fibre metal laminates. *Materials & Design* **98**, 98-107 (2016).

20. F. Sarasini, J. Tirillo, S. D'Altilia, T. Valente, C. Santulli, F. Touchard, L. Chocinski-Arnault, D. Mellier, L. Lampani and P. Gaudenzi, Damage tolerance of carbon/flax hybrid composites subjected to low velocity impact. *Composite Part B* **91**, 144-153 (2016).

21. E. Nisini, C. Santulli, A. Liverani, Mechanical and impact characterization of hybrid composite laminates with carbon, basalt and flax fibres. *Composite Part B* **127**, 92-99 (2017).

22. H.N. Dhakal, Z.Y. Zhang, N. Bennett and P.N.B. Reis, Low-velocity impact response of non-woven hemp fibre reinforced unsaturated polyester composites: Influence of impactor geometry and impact velocity. *Journal of Composite Structures* **4**, 2756-2763 (2012).

23. G. Caprino, L. Carrino, M. Durante, A. Langella and V. Lopresto, Low impact 475ehavior of hemp fibre reinforced epoxy composites. *Journal of Composite Structures* **133**, 892-901 (2015).

24. C. Scarponi, F. Sarasini, J. Tirillo and L. Lampani, Low-velocity impact behaviour of hemp fibre reinforced bio-based epoxy laminates. *Composite Part B* **91**, 162-168 (2016).

25. S. Liang, L. Guillaumat and P. B. Gning, Impact behavior of flax/epoxy composite plates. *International Journal of Impact Engineering* **80**, 56-64 (2015).

26. R. Hill, A theory of the yielding and plastic flow of anisotropic metals. *Proceedings of the Royal Society A* **193**, 281-297 (1948).

27. M. Wali, H. Chouchene, L. Ben Said and F. Dammak, One-equation integration algorithm of a generalized quadratic yield function with Chaboche non-linear isotropic/kinematic hardening. *International Journal of Mechanical Sciences* **92**, 223-232 (2015).

28. J. Mars, S. Koubaa, M. Wali and F. Dammak, Numerical analysis of geometrically non-linear behavior of functionally graded shells. *Latin American Journal of Solids and Structures* **14**(11), 1952-1978 (2017).

29. J. Mars, L. Ben Said, M. Wali and F. Dammak, Elasto-plastic modelling of low-velocity impact on functionally graded circular plates. *International Journal of Applied Mechanics* **10**(4) (2018).

30. M. Wali, R. Autay, J. Mars and F. Dammak, A simple integration algorithm for a non-associated anisotropic plasticity model for sheet metal forming. *International Journal of Mechanical Engineering* **10**(11), 1952-1978 (2017).
REFERENCES

1. H.L. Bos, The potential of flax fibres as reinforcement for composite materials (Ph.D. Thesis). Eindhoven: Technische Universiteit Eindhoven (2004).

2. C. Bailey, Analysis of the flax fibres tensile behavior and analysis of the tensile stiffness increase. *Composite Part A* **33**, 939-948 (2002).

3. C. Gourier, A. Le Duigou, A. Bourmaud and C. Bailey, Mechanical analysis of elementary flax fibre tensile properties after different thermal cycles. *Composite Part A* **64**, 159-166 (2014).

4. M. Fehri, R.R. Ragueh, A. Vivet, F. Dammak and M. Haddar, Improvement of natural fiber composite materials by carbon fibers. *Journal of Renewable Materials* **5**(1), 38-47 (2017).

5. J. Andersons, J. Modniks and E. Sparnins, Modeling the nonlinear deformation of flax-fiber-reinforced polymer matrix laminates in active loading. *Journal of Reinforced Plastics and Composites* **34**(3), 248-256 (2015).

6. A. Rubio-López, T. Hoang and C. Santiuste, Constitutive model to predict the viscoplastic behaviour of natural fibres based composites. *Journal of Composite Structures* **155**, 8-18 (2016).

7. C. Poilane, Z.E. Cherif, F. Richard, A. Vivet, B. Ben Doudou and J. Chen, Polymer reinforced by flax fibres as a viscoelastoplastic material. *Journal of Composite Structures* **112**, 100-112 (2014).

8. Z.E. Cherif, C. Poilâne, A. Vivet, B. Ben Doudou and J. Chen. About optimal architecture of plant fibre textile composite for mechanical and sorption properties. *Journal of Composite Structures* **140**, 240-251 (2016).

9. M.J. John and S. Thomas, Biofibres and biocomposites. *Carbohydr Polym* **71**, 343-364 (2008).

10. A.K. Bledzki and A. Jaszkiewicz, Mechanical performance of biocomposites based on PLA and PHBV reinforced with natural fibres-a comparative study to PP. *Composites Science and Technology* **70**, 1687-1696 (2010).

11. H. Ku, Y.M. Cheng, C. Snook and D. Baddeley, Drop weight impact test fracture of vinyl ester composites: micrographs of pilot study. *Journal of Composite Materials* **39**, 1607-1620 (2005).

12. C. Santiuste, S. Sanchez-Saez and E. Barbero, Residual flexural strength after low velocity impact in glass/polyester composite beams. *Journal of Composite Structures* **92**, 25-30 (2010).

13. C. Santiuste, S. Sanchez-Saez and E. Barbero, A comparison of progressive-failure criteria in the prediction of the dynamic bending failure of composite laminated beams. *Journal of Composite Structures* **92**, 2406-2414 (2010).

14. I. Ivañez, C. Santiuste and S. Sanchez-Saez, FEM analysis of dynamic flexural behavior of composite sandwich beams with foam core. *Journal of Composite Structures* **92**, 2285-2291 (2010).

15. A.R. Setoodeh, P. Malekzadeh and K. Nikbin, Low velocity impact analysis of laminated composite plates using a 3D elasticity based layerwise FEM. *Materials & Design* **30**, 3795-3801 (2009).

16. S.H. Song, Y.S. Byun, T.W. Ku, W.J. Song, J. Kim and B.S. Kang, Experimental and numerical investigation on impact performance of carbon reinforced aluminum laminates. *Journal of Materials Science & Technology* **26**, 327-332 (2010).

17. S.M.R. Khalili, M. Soroush, A. Davar and O. Rahmani, Finite element modeling of low-velocity impact on laminated composite plates and cylindrical shells. *Journal of Composite Structures* **93**, 1363-1375 (2011).

18. D. Zhang, Y. Sun, L. Chen and N. Pan, A comparative study on low-velocity impact response of fabric composite laminates. *Materials & Design* **50**, 750-756 (2013).

19. L. Ferrante, F. Sarasini, J. Tirillò, L. Lampani, T. Valente and P. Gaudenzi, Low velocity impact response of basalt-aluminium fibre metal laminates. *Materials & Design* **98**, 98-107 (2016).

20. F. Sarasini, J. Tirillo, S. D'Altilia, T. Valente, C. Santulli, F. Touchard, L. Chocinski-Arnault, D. Mellier, L. Lampani and P. Gaudenzi, Damage tolerance of carbon/flax hybrid composites subjected to low velocity impact. *Composite Part B* **91**, 144-153 (2016).

21. E. Nisini, C. Santulli, A. Liverani, Mechanical and impact characterization of hybrid composite laminates with carbon, basalt and flax fibres. *Composite Part B* **127**, 92-99 (2017).

22. H.N. Dhakal, Z.Y. Zhang, N. Bennett and P.N.B. Reis, Low-velocity impact response of non-woven hemp fibre reinforced unsaturated polyester composites: Influence of impactor geometry and impact velocity. *Journal of Composite Structures* **4**, 2756-2763 (2012).

23. G. Caprino, L. Carrino, M. Durante, A. Langella and V. Lopresto, Low impact 475ehavior of hemp fibre reinforced epoxy composites. *Journal of Composite Structures* **133**, 892-901 (2015).

24. C. Scarponi, F. Sarasini, J. Tirillo and L. Lampani, Low-velocity impact behaviour of hemp fibre reinforced bio-based epoxy laminates. *Composite Part B* **91**, 162-168 (2016).

25. S. Liang, L. Guillaumat and P. B. Gning, Impact behavior of flax/epoxy composite plates. *International Journal of Impact Engineering* **80**, 56-64 (2015).

26. R. Hill, A theory of the yielding and plastic flow of anisotropic metals. *Proceedings of the Royal Society A* **193**, 281-297 (1948).

27. M. Wali, H. Chouchene, L. Ben Said and F. Dammak, One-equation integration algorithm of a generalized quadratic yield function with Chaboche non-linear isotropic/kinematic hardening. *International Journal of Mechanical Sciences* **92**, 223-232 (2015).

28. J. Mars, S. Koubaa, M. Wali and F. Dammak, Numerical analysis of geometrically non-linear behavior of functionally graded shells. *Latin American Journal of Solids and Structures* **14**(11), 1952-1978 (2017).

29. J. Mars, L. Ben Said, M. Wali and F. Dammak, Elasto-plastic modelling of low-velocity impact on functionally graded circular plates. *International Journal of Applied Mechanics* **10**(4) (2018).

30. M. Wali, R. Autay, J. Mars and F. Dammak, A simple integration algorithm for a non-associated anisotropic plasticity model for sheet metal forming. *International Journal of Plasticity* **100**, 1-12 (2000).
*Journal for Numerical Methods in Engineering* **107**(3), 183-204 (2016).

31. J.C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. *Computer Methods in Applied Mechanics and Engineering* **99**, 161-112 (1992).

32. S. Koubaa, J. Mars, M. Wali and F. Dammak, Numerical study of anisotropic behavior of Aluminum alloy subjected to dynamic perforation. *International Journal of Impact Engineering* **101**, 105-114 (2017).

33. L. Belhassen, S. Koubaa, M. Wali and F. Dammak, Numerical prediction of springback and ductile damage in rubber-pad forming process of aluminum sheet metal. *International Journal of Mechanical Sciences* **117**, 218-226 (2016).

34. L. Ben Said, J. Mars, M. Wali and F. Dammak, Numerical prediction of the ductile damage in single point incremental forming process. *International Journal of Mechanical Sciences* **131-132**, 546-558 (2017).

35. R. Autay, S. Koubaa, M. Wali and F. Dammak, Numerical implementation of coupled anisotropic plasticity-ductile damage in sheet metal forming process. *Journal of Mechanics* **34**(4), 417-430 (2017).

36. J. Berthe, M. Brieu and E. Deletombe, Thermo-viscoelastic modeling of organic matrix composite behavior-Application to T700GC/M21. *Mechanics of Materials* **81**, 18-24 (2015).

37. S. Treutenaere, F. Lauro, B. Bennani, G. Haugou, T. Matsumoto and E. Mottola, Constitutive modeling of the strain-rate dependency of fabric reinforced polymers. *International Journal of Impact Engineering* **108**, 361-369 (2017).

38. P. Ladeveze, O. Allix, J. F. Deu and D. Léveque, A mesomodel for localization and damage computation in laminates. *Computer Methods in Applied Mechanics and Engineering* **183**, 105-122 (2000).

39. J. Mars, M. Wali, A. Jarraya, F. Dammak and A. Dhiab, Finite element implementation of an orthotropic plasticity model for sheet metal in low velocity impact simulations. *Thin Wall Structures* **89**, 93-100 (2015).

40. J. Mars, E. Chebbi, M. Wali and F. Dammak, Numerical and experimental investigations of low velocity impact on glass fiber-reinforced polyamide. *Composite Part B: Engineering* **146**, 116-123 (2018).

41. M. Wali, M. Abdennadher, T. Fakhfakh and M. Haddar, Dynamic analysis of an elasto-plastic sandwich subjected to low velocity impact. *Multidiscipline Modeling in Materials and Structures* **7**(2), 184-206 (2011).
*Journal for Numerical Methods in Engineering* **107**(3), 183-204 (2016).

31. J.C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. *Computer Methods in Applied Mechanics and Engineering* **99**, 161-112 (1992).

32. S. Koubaa, J. Mars, M. Wali and F. Dammak, Numerical study of anisotropic behavior of Aluminum alloy subjected to dynamic perforation. *International Journal of Impact Engineering* **101**, 105-114 (2017).

33. L. Belhassen, S. Koubaa, M. Wali and F. Dammak, Numerical prediction of springback and ductile damage in rubber-pad forming process of aluminum sheet metal. *International Journal of Mechanical Sciences* **117**, 218-226 (2016).

34. L. Ben Said, J. Mars, M. Wali and F. Dammak, Numerical prediction of the ductile damage in single point incremental forming process. *International Journal of Mechanical Sciences* **131-132**, 546-558 (2017).

35. R. Autay, S. Koubaa, M. Wali and F. Dammak, Numerical implementation of coupled anisotropic plasticity-ductile damage in sheet metal forming process. *Journal of Mechanics* **34**(4), 417-430 (2017).

36. J. Berthe, M. Brieu and E. Deletombe, Thermo-viscoelastic modeling of organic matrix composite behavior-Application to T700GC/M21. *Mechanics of Materials* **81**, 18-24 (2015).

37. S. Treutenaere, F. Lauro, B. Bennani, G. Haugou, T. Matsumoto and E. Mottola, Constitutive modeling of the strain-rate dependency of fabric reinforced polymers. *International Journal of Impact Engineering* **108**, 361-369 (2017).

38. P. Ladeveze, O. Allix, J. F. Deu and D. Léveque, A mesomodel for localization and damage computation in laminates. *Computer Methods in Applied Mechanics and Engineering* **183**, 105-122 (2000).

39. J. Mars, M. Wali, A. Jarraya, F. Dammak and A. Dhiab, Finite element implementation of an orthotropic plasticity model for sheet metal in low velocity impact simulations. *Thin Wall Structures* **89**, 93-100 (2015).

40. J. Mars, E. Chebbi, M. Wali and F. Dammak, Numerical and experimental investigations of low velocity impact on glass fiber-reinforced polyamide. *Composite Part B: Engineering* **146**, 116-123 (2018).

41. M. Wali, M. Abdennadher, T. Fakhfakh and M. Haddar, Dynamic analysis of an elasto-plastic sandwich subjected to low velocity impact. *Multidiscipline Modeling in Materials and Structures* **7**(2), 184-206 (2011).