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

HAL
open science

A rate dependent constitutive model for
glass-fiber/epoxy-matrix woven fabrics

Steven Marguet, Patrick Rozycki, Laurent Gornet

► To cite this version:

Steven Marguet, Patrick Rozycki, Laurent Gornet. A rate dependent constitutive model for glass-fiber/epoxy-matrix woven fabrics. CMC-Computers, Materials & Continua, 2006, 4 (3), pp.119-136. 10.3970/cmc.2006.004.119. hal-02151389

HAL Id: hal-02151389

https://hal.science/hal-02151389

Submitted on 8 Jun 2019

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

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

Public Domain
# A Rate Dependent Constitutive Model for Glass-Fibre / Epoxy-Matrix Woven Fabrics

S. Marguet<fnref n="1" />, P. Rozycki<fnref n="1" /> and L. Gornet<fnref n="1" />

**Abstract:** This paper deals with the modelling until rupture of composite structures made of glass-fibre / epoxy-resin woven fabrics submitted to dynamic loadings. The model is built at the mesoscale of the elementary ply. It takes into account the slightly nonlinear brittle behaviour of the fibres under tensile sollicitations, their nonlinear behaviour in compression as well as the strongly non linear and irreversible behaviour of the ply in shear. Strain rate effects are also introduced and special attention is paid to the objectivity of the model in the context of finite element calculation. Therefore the choice of a delayed damage mesomodel coupled with viscoplasticity is made. In order to identify the values of the parameters of the model, an optimization procedure based on a gradient-free direct search method has been developed. As a logical proceeding to this study, the model ability to avoid strain localization and mesh dependence is then checked on simple uniaxial examples. The last part of this paper is devoted to structural calculation. The results of the simulations of both the impact on composite plate and the crushing of thin-walled tube demonstrate the capability of the model to reproduce observed physical phenomena.

**keyword:** Composite woven fabrics, delayed damage mesomodel, viscoplasticity, identification, direct search method, dynamic loadings.

## 1 Introduction

Nowadays composite materials take a growing importance in the design and the manufacture of multihull oceanic sailing ships. Their ability to provide significant stiffness combined with limited weight perfectly matches the more and more constraining needs of competition. Designing the ships with reduced safety margins is also a good way to improve the performances measured in term

of the "velocity / wind power" ratio. Still, during great oceanic races, structural failures caused by extreme climatic conditions lead lots of ships to give up the race. The knowledge of the material behaviour under dynamic loadings such as slamming or wave impacts then appears to be determining for the design of safe structures. Fibre reinforced polymer resin woven fabrics composite materials have been studied for many years. Their linear orthotropic elastic behaviour is now well known as can be seen in [Herakovich (1998)]. The use of homogenization procedures to find the equivalent mechanical properties of stacking sequences of these materials has allowed the development of huge numerical models dedicated to the design of the sailing ships [Gornet, Marguet, and Marckmann (2006)]. However, in a context of competition, the knowledge of the elastic properties of the constitutive material is not enough for an optimized design of the structures. Nonlinear effects which correspond to the local deteriorations of the composite materials and eventually lead to the failure of the whole ship have to be dealt with.

Since the end of the 80's, an important work was performed on the understanding and modelling of the damage effects. Four principal modes of degradation of the elementary ply were isolated:

* the tensile fragile rupture of fibres in tension;

* the local buckling of the fibres in compression (responsible of the observed nonlinearity);

* the microcracking of the resin under shear or transverse loadings (which explains the nonlinearity and the irreversible strains);

* and the decohesion between the fibres and the matrix.

To these modes can be added the delamination between the different staked layers of the composite which also plays an important role in the failure of macroscopic

<footnote n="1">¹ steven.marguet@ec-nantes.fr, Ecole Centrale de Nantes, GeM, UMR CNRS 6183, 1 rue de la Noë, BP 92101, 44321 Nantes Cedex 03, France.</footnote>
structures. Based upon these observations, several models written at the homogenized mesoscale of the elementary ply have been developed and especially by Ladevèze and his co-workers [Ladevèze and LeDantec (1992), Ladevèze, Allix, Deü, and Lévêque (2000), Ladevèze, Allix, Gornet, Lévêque, and Perret (1998)]. In a first step, damage mesomodels have demonstrated their ability to take into account the loss of stiffness observed in experiments. The adjonction of plastic models with isotropic hardening has also enabled to represent the irreversible strains. Yet, in the context of finite element calculation, local softening damage models are subjected to the strain localization phenomenon which results in a state of damage dependent from the elements. To circumvent this issue, delayed damage mesomodels were developed by [Ladevèze, Allix, Deü, and Lévêque (2000)]. The aim is to introduce an internal length linked with a maximal rate of damage. As a consequence, the damage can not immediately reach its critical value, and so, information can be passed through the mesh. Numerous studies have then been performed in order to check the ability of delayed damage mesomodels to avoid strain localization [Suffis, Lubrecht, and Combescure (2003)]. By the way, the identification of the values of the parameters and the validity of such models on large strain rate ranges still remain a matter. More recently, a model resting on classical damage and plastic theories was developed by [Rozycki (2000)] to deal with the strain rate dependency of the glass-fibre / epoxy-resin composite materials. Satisfying results were obtained but the mesh independence was not ensured. The object of this study is to propose a local mesomodel dedicated to glass-fibre / epoxy-resin woven fabrics submitted to high strain rate loadings. An optimization procedure based on a gradient-free direct search method (the Pattern Search Algorithm) is also implemented in order to perform the identification of the parameters of the model. To check the behaviour of the model in finite element calculations, the example of uniaxial tensile tests on bars was taken. The ability of the model to reproduce the failure of composite structures was eventually showed thanks to the simulation of the impact on a woven fabric composite plate, and of the crushing of a thin-walled composite tube.

## 2 Modelling of the fabric

The present model is written at the mesoscale of the layer which brings a good compromise between the scales of

both the constitutive materials (glass fibres and epoxy resin) and the structure. Moreover, and thanks to the thin aspect of most of the composite structures, the assumption of plane stress is done. At last, the effects of temperature are neglected and only small strains are considered. From now on, the subscripts 1 and 2 stand for the orthogonal fibres directions while 3 stands for the direction normal to the woven fabric ply.

### 2.1 Thermodynamic potential and state laws

To deal with both the reversible and the irreversible effects, the classical split of the total train $\varepsilon^t$ is assumed:

$$ \varepsilon^t = \varepsilon^e + \varepsilon^{vp} \quad (1) $$

where $\varepsilon^e$ stands for the elastic strain which is an internal variable and $\varepsilon^{vp}$ stands for the viscoplastic strain. As viscoplasticity only occurs in shear, where fibres cannot prevent for viscoplastic flow, the different strain vectors take the form $\varepsilon^t = [\varepsilon_{11}^t \ \varepsilon_{22}^t \ \varepsilon_{12}^t]^T$; $\varepsilon^e = [\varepsilon_{11}^e \ \varepsilon_{22}^e \ \varepsilon_{12}^e]^T$ and, for irreversible strains, $\varepsilon^{vp} = [0 \ 0 \ \varepsilon_{12}^{vp}]^T$.

To take into account the irreversible strains, the loss of rigidity, the different tensile and compressive behaviour in fibre direction and the strain rate dependency, the thermodynamic potential is derived from the Helmholtz's volumic free energy and written as follows:

$$ \begin{aligned} 2 \cdot \psi = 2 \cdot \psi(\varepsilon^e; \underline{d}; p; (\dot{\varepsilon}^e; t)) \\ = & \left( (\varepsilon_{11}^e)^2 \cdot (1 - d_{11}) \cdot \Lambda_{11}^0 \cdot V_{11} (\dot{\varepsilon}_{11}^e) \right) \cdot P^+ (\sigma_{11}^{t-\Delta t}) \\ & + \left( (\varepsilon_{11}^e)^2 \cdot (1 - d_{11}^c) \cdot \Lambda_{11}^0 \cdot V_{11} (\dot{\varepsilon}_{11}^e) \right) \cdot P^- (\sigma_{11}^{t-\Delta t}) \\ & + \varepsilon_{11}^e \cdot \Lambda_{12}^0 \cdot V_{12} (\dot{\varepsilon}_{12}^e; \dot{\varepsilon}_{22}^e) \cdot \varepsilon_{22}^e \\ & + \left( (\varepsilon_{22}^e)^2 \cdot (1 - d_{22}) \cdot \Lambda_{22}^0 \cdot V_{22} (\dot{\varepsilon}_{22}^e) \right) \cdot P^+ (\sigma_{22}^{t-\Delta t}) \\ & + \left( (\varepsilon_{22}^e)^2 \cdot (1 - d_{22}^c) \cdot \Lambda_{22}^0 \cdot V_{22} (\dot{\varepsilon}_{22}^e) \right) \cdot P^- (\sigma_{22}^{t-\Delta t}) \\ & + \varepsilon_{22}^e \cdot \Lambda_{21}^0 \cdot V_{21} (\dot{\varepsilon}_{11}^e; \dot{\varepsilon}_{22}^e) \cdot \varepsilon_{11}^e \\ & + (2 \cdot \varepsilon_{12}^e)^2 \cdot (1 - d_{12}) \cdot \Lambda_{1212}^0 \cdot V_{1212} (\dot{\varepsilon}_{12}^e) + h(p) \end{aligned} \quad (2) $$

where:

* $\dot{\varepsilon}^e$ must be considered as a parameter, it allows the elastic modulus to increase with the strain rate as observed in experiments;

* $\Lambda_{ii}^0$ are the coefficients of the stiffness matrix associated with the undamage material which is supposed to be the initial one: $\Lambda_{11}^0 =$
$\frac{E_{11}}{1-v_{12}v_{21}}$ ; $\Lambda_{22}^{0} = \frac{E_{22}}{1-v_{12}v_{21}}$ ; $\Lambda_{12}^{0} = \frac{v_{21} \cdot E_{11}}{1-v_{12}v_{21}} = \Lambda_{21}^{0} = \frac{v_{12} \cdot E_{22}}{1-v_{12}v_{21}}$ ; $\Lambda_{1212}^{0} = G_{12}$ with $E_{11}, E_{22}$ and $G_{12}$ the elastic moduli and $v_{12}, v_{21}$ the Poisson's moduli.

* $d_{11}$ and $d_{22}$ are damage variables that drive the brittle rupture of the fibres under tensile sollicitations;

* $d_{11}^c$ and $d_{22}^c$ are pseudo damage variables that stand for the nonlinearity under compressive state due to the microbuckling of the fibres, at time $t$, they are explicitly expressed as a function of the stresses: $d_{ii}^c = \rho_{ii} \cdot \sigma_{ii}^{j-\Delta t}$ for $ii = \{11; 22\}$ and with $\rho_{ii}$ material parameters;

* $d_{12}$ is a damage variable associated with the loss of stiffness in shear due to the matrix microcracking;

* and $h(p)$ is the volumic energy associated with a nonlinear isotropic hardening ($p$ is the cumulated plastic strain).

To complete the description of the state potential, the following functions have been introduced:

* $P^+(\sigma^{t-\Delta t}) = 1$ if $\sigma^{t-\Delta t} \ge 0$, 0 otherwise (this is an explicit pragmatic way to check for the tensile or compressive state of the material from the state of stress at one step time before the actual time $t$);

* $P^-(\sigma^{t-\Delta t}) = 1$ if $\sigma^{t-\Delta t} < 0$, 0 otherwise;

* $V_{ii}(\dot{\epsilon}_{ii}^e) = (1 + v_{ii} \cdot \ln(\frac{|\dot{\epsilon}_{ii}^e|}{\dot{\epsilon}_{ref}^e}))$ for $ii = \{11; 22; 1212\}$ and for $\dot{\epsilon}_{ii}^e \in [\dot{\epsilon}_{min}^e; \dot{\epsilon}_{max}^e]$ are functions that drive the rate dependency of the stiffness of the material;

* and $V_{12}(\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) = V_{21}(\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) = (1 + \ln(\frac{|\dot{\epsilon}_{11}^e|^{v_{11}} + |\dot{\epsilon}_{22}^e|^{v_{22}}}{\dot{\epsilon}_{ref}^{v_{11}+v_{22}}}))^{\frac{1}{2}}$ for $(\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \in [\dot{\epsilon}_{min}^e; \dot{\epsilon}_{max}^e]^2$

do the same for the nondiagonal terms of the stiffness matrix and allows to conserve its symmetry.

Here it has to be noticed that the introduction of all the preceeding functions preserves the symmetry and the positive definite property of the stiffness matrix.

The state laws can now be derived from equation (2) leading to the expressions of $\underline{\sigma}$, the Cauchy's stress vector dual of the elastic strain, $\underline{Y}$ the thermodynamic forces

associated with the damage state and $R$ the stress linked to the isotropic hardening:

$$ \underline{\sigma} = \frac{\partial \psi}{\partial \underline{\epsilon}^e}; \underline{Y} = -\frac{\partial \psi}{\partial \underline{d}}; \underline{R} = \frac{\partial \psi}{\partial p} \quad (3) $$

Taken the volumic energy due to isotropic hardening under the form $h(p) = \frac{Q}{\beta+1} \cdot p^{\beta+1}$, the preceeding expressions once developed, the state laws become:

$$ \begin{aligned} \sigma_{11} = & ((1-d_{11}) \cdot \Lambda_{11}^0 \cdot V_{11} (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e) \cdot P^+ (\sigma_{11}^{t-\Delta t}) \\ & + ((1-d_{11}^c) \cdot \Lambda_{11}^0 \cdot V_{11} (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e) \cdot P^- (\sigma_{11}^{t-\Delta t}) \\ & + \Lambda_{12}^0 \cdot V_{12} (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e \end{aligned} $$

$$ \begin{aligned} \sigma_{22} = & ((1-d_{22}) \cdot \Lambda_{22}^0 \cdot V_{22} (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e) \cdot P^+ (\sigma_{22}^{t-\Delta t}) \\ & + ((1-d_{22}^c) \cdot \Lambda_{22}^0 \cdot V_{22} (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e) \cdot P^- (\sigma_{22}^{t-\Delta t}) \\ & + \Lambda_{21}^0 \cdot V_{21} (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{11}^e \end{aligned} $$

$$ \sigma_{12} = (1 - d_{12}) \cdot \Lambda_{1212}^0 \cdot V_{1212} (\dot{\epsilon}_{12}^e) \cdot (2 \cdot \epsilon_{12}^e) $$

$$ Y_{11} = \frac{1}{2} \cdot \Lambda_{11}^0 \cdot V_{11} (\dot{\epsilon}_{11}^e) \cdot (\epsilon_{11}^e)^2 $$

$$ Y_{22} = \frac{1}{2} \cdot \Lambda_{22}^0 \cdot V_{22} (\dot{\epsilon}_{22}^e) \cdot (\epsilon_{22}^e)^2 $$

$$ Y_{12} = 2 \cdot \Lambda_{1212}^0 \cdot V_{1212} (\dot{\epsilon}_{12}^e) \cdot (\epsilon_{12}^t - \epsilon_{12}^{vp})^2 $$

$$ R = Q \cdot p^{\beta} \quad (4) $$

## 2.2 Dissipation potential and evolution laws

To obtain the evolution laws, the existence of a dissipation potential that is a function of the stress variables is assumed [Lemaître and Chaboche (1990)]. The following expression is adopted in order to deal with both the rate dependent irreversible effects and the degradations:

$$ \phi(\underline{\sigma}; \underline{Y}; \underline{R}; (\underline{d})) = \phi^{vp}(\underline{\sigma}; \underline{Y}; \underline{R}; (\underline{d})) + \phi^{d}(\underline{\sigma}; \underline{Y}; \underline{R}; (\underline{d})) \quad (5) $$

where $\phi^{vp}$ stands for the dissipated energy per unit volume due to viscoplasticity and $\phi^d$ is the volumic energy dissipated during the microcracking of fibres and resin.
Introducing the Perzyna's viscoplastic theory [Perzyna (1966)], and remembering that irreversible strains only occur in shear, the viscoplastic dissipation potential is:

$$ \Phi^{vp}(\sigma_{12}; R; (d_{12})) = \frac{K}{m+1} \left\langle \frac{f(\sigma_{12}; R; (d_{12}))}{K} \right\rangle_{+}^{m+1} \quad (6) $$

with $f$, the yield surface, $K$ and $m$ parameters that drive the rate sensitivity of the material, and $d_{12}$ the damage variable associated to the microcracking of the resin in shear. The choice of such a potential is motivated by two different aspects:

* the dependence of the woven fabric to the loading rate (material aspect that should come from the viscoelastic nature of the epoxy resin);

* and the viscoplastic regularization effects on strain localization phenomenon (numerical aspect).

The coupling between viscoplasticity and damage is realized through the introduction of the effective stress in the expression of the yield surface $f$:

$$ f(\sigma_{12}; R; (d_{12})) = f(\tilde{\sigma}_{12}; R) = |\tilde{\sigma}_{12}| - R - \sigma^v \quad (7) $$

with $\tilde{\sigma}_{12} = \frac{\sigma_{12}}{1-d_{12}}$ the effective shear stress which stands for the stress that should be applied to the undamage material to obtain the same state of deformation as the one obtained with $\sigma_{12}$ on the damaged material.

From equation (6) it is now possible to determine the evolution laws of both the inelastic strain $\epsilon_{12}^{vp}$ and the cumulated plastic strain $p$:

$$ \begin{aligned} \dot{\epsilon}_{12}^{vp} &= \frac{\partial \phi}{\partial \sigma_{12}} = \dot{\lambda}^{vp} \frac{\partial f}{\partial \sigma_{12}} \\ \dot{p} &= -\frac{\partial \phi}{\partial R} = -\dot{\lambda}^{vp} \frac{\partial f}{\partial R} \end{aligned} \quad (8) $$

where:

$$ \dot{\lambda}^{vp} = \left\langle \frac{f}{K} \right\rangle_{+}^{m} \quad (9) $$

corresponds to a Lagrange's multiplier that drives the viscoplastic flow and is named viscoplastic multiplier. The function $\langle x \rangle_+$ stands for the Heaviside's function: $\langle x \rangle_+ = x$ if $x > 0$ and $0$ elsewhere. If $f < 0$ no viscoplastic flow can occur and then $\dot{\epsilon}_{12}^{vp}$ and $\dot{p}$ are both equal to

zero. If $f > 0$ then viscoplasticity appears and the material flows in the direction normal to the elastic surface.

To completely define the dissipation potential $\phi$, the choice of the damage dissipation potential $\phi^d$ has to be carried out. Assumption is made that a dissipation potential, depending only on the damage state, exists. The following evolution laws for the damage variables are then obtained:

$$ \dot{d}_{ii} = \frac{1}{\tau_{ii}} \cdot \left(1 - e^{-\alpha_{ii} \cdot \langle z_{ii} \rangle_{+}}\right) \quad (10) $$

if $d < d_{ii}^c$; 0 elsewhere, with $ii = \{11; 22; 12\}$ and where $\alpha_{ii}$ and $\tau_{ii}$ are material parameters and $z_{ii}$ stand for the functions that drive the appearance of the damage and its rate.

In fibre direction, $z_{ii} = \frac{\sqrt{Y_{ii}} - \sqrt{Y_{ii}^o}}{\sqrt{Y_{ii}^c} - \sqrt{Y_{ii}^o}} - d_{ii}$ ($ii = 11, 22$) whereas in shear direction

$z_{12} = \frac{\ln(Y_{12}) - \ln(Y_{12}^o)}{\ln(Y_{12}^c) - \ln(Y_{12}^o)} - d_{12}$ with $Y_{ii}^o$ the yield thermodynamic force that initiates the damage and $Y_{ii}^c$ a critical thermodynamic force.

The evolution equations of the damage variables $d$ define a damage mesomodel as was introduced by Ladevèze and his coworkers [Ladevèze, Allix, Deü, and Lévêque (2000)]. It is particularly interesting for its ability to avoid the strain localization phenomenon which results in mesh dependence of the results (see [Suffis, Lubrecht, and Combescure (2003), Chung, Choi, and Kim (2003)] for more details). The material parameter $\tau_{ii}$ represents the inverse of the maximum damage rate and is linked to an internal length supposed to be the thickness of the layer for the mesomodel. As a consequence, the damage has to be constant through the out-of-plan direction of a layer of the composite ply.

A last key point in the construction of the model is to check the positivity of the intrinsic dissipation, coming from the second law of thermodynamics, and defined by the equation:

$$ D = \underline{\sigma}^T \dot{\underline{\epsilon}}^t - \dot{\psi} = \underline{\sigma}^T \dot{\underline{\epsilon}}^{vp} - \underline{Y}^T \dot{\underline{d}} - R \dot{p} \quad (11) $$

As no simple expression can be obtained at this point, particular attention should be paid to the verification of the second thermodynamic law.

## 2.3 Implementation and numerical aspects

This section deals with the following problem: "given the state of the material on the whole model at
time $t_n$, known the step time $\Delta t_{n+1} = t_{n+1} - t_n$ and the prescribed conditions $\underline{\underline{\epsilon}}_{n+1}^t$, the state of the material at time $t_{n+1}$ is to be found".

The integration is performed in two major steps which are:

1. integration of the delayed damage model in fibre directions;

2. integration of the coupled viscoplastic/delayed damage mesomodel in shear thanks to a modified Radial Return Mapping algorithm [Simo and Hugues (1998)].

For the behaviour in fibre direction, no special difficulty occurs so that the integration method is not detailed. The damage state at time $t_{n+1}$ is simply computed implicitly by the mean of a Newton Raphson's procedure considering whether the type of sollicitation at time $t_n$ is tension or compression. At the end of this first part, the damage variables $d_{11}^t$ and $d_{22}^t$ are updated and are used to check for the failure of the fibres in tension. Under compressive state, a simple criterion based on a critical stress drives the rupture.

Next step deals with the integration of the coupled viscoplastic/delayed damage mesomodel. The system of equations to solve comes from the implicit backward Euler's scheme applied to the variables from strain space ($d, \dot{d}, \epsilon^{vp}$ and $\dot{\epsilon}^{vp}$ stand for $d_{12}, \dot{d}_{12}, \epsilon_{12}^{vp}$ and $\dot{\epsilon}_{12}^{vp}$ respectively):

$$d_{n+1} = d_n + \Delta t_{n+1} \dot{d}_{n+1}$$

$$\epsilon_{n+1}^{vp} = \epsilon_{n+1}^{vp} + \Delta t_{n+1} \dot{\epsilon}_{n+1}^{vp}$$

$$p_{n+1} = p_n + \Delta t_{n+1} \dot{p}_{n+1} \quad (12)$$

where the evolution laws give the last right hand term of the equations. To condense these expressions, it is introduced the vector of unknowns $\underline{\underline{\Theta}} = (d; \epsilon^{vp}; p)^T$ which leads to:

$$\underline{\underline{\Theta}}_{n+1} = \underline{\underline{\Theta}}_n + \Delta t_{n+1} \dot{\underline{\underline{\Theta}}}_{n+1} \quad (13)$$

As the expression of $\dot{\underline{\underline{\Theta}}}$ is a non linear expression of $\underline{\underline{\Theta}}$, the solution of the system is not trivial. In order to solve it, an incremental procedure of Newton Raphson is once again employed. At the increment $k$ of the iterative procedure a solution to the problem at time $t_{n+1}$, noted $\underline{\underline{\Theta}}_{n+1}^k$,

is supposed to be known (for the first iteration $k=0$, the approximation $\underline{\underline{\Theta}}_{n+1}^0 = \underline{\underline{\Theta}}_n$ is used). This hypothesis leads to the formation of a residual vector:

$$\underline{\underline{Res}}_{n+1}^k = \underline{\underline{\Theta}}_{n+1}^k - \underline{\underline{\Theta}}_n - \Delta t_{n+1} \dot{\underline{\underline{\Theta}}}_{n+1}^k \quad (14)$$

If the norm of the residual is greater than a given tolerance, the solution is not acceptable, and an iteration has to be performed. Linearizing the equation (14) conducts to the following expression of the increment in strain variables:

$$\delta \underline{\underline{\Theta}}_{n+1}^k = - \left( \left[ \frac{\partial \underline{\underline{Res}}}{\partial \underline{\underline{\Theta}}} \right]_{n+1}^k \right)^{-1} \underline{\underline{Res}}_{n+1}^k \quad (15)$$

with

$$\left[ \frac{\partial \underline{\underline{Res}}}{\partial \underline{\underline{\Theta}}} \right]_{n+1}^k = \begin{bmatrix} 1 - \Delta t_{n+1} \frac{\partial \dot{d}}{\partial d} & - \Delta t_{n+1} \frac{\partial \dot{d}}{\partial \epsilon^{vp}} & - \Delta t_{n+1} \frac{\partial \dot{d}}{\partial p} \\ - \Delta t_{n+1} \frac{\partial \dot{\epsilon}^{vp}}{\partial d} & 1 - \Delta t_{n+1} \frac{\partial \dot{\epsilon}^{vp}}{\partial \epsilon^{vp}} & - \Delta t_{n+1} \frac{\partial \dot{\epsilon}^{vp}}{\partial p} \\ - \Delta t_{n+1} \frac{\partial \dot{p}}{\partial d} & - \Delta t_{n+1} \frac{\partial \dot{p}}{\partial \epsilon^{vp}} & 1 - \Delta t_{n+1} \frac{\partial \dot{p}}{\partial p} \end{bmatrix}_{n+1}^k \quad (16)$$

which clearly shows the coupling between viscoplasticity and delayed damage. From the preceding two equations, it is then possible to infer a new estimation of the state of the material at time $t_{n+1}$ with:

$$\underline{\underline{\Theta}}_{n+1}^k = \underline{\underline{\Theta}}_{n+1}^k + \delta \underline{\underline{\Theta}}_{n+1}^k \quad \text{and} \quad k = k+1 \quad (17)$$

which leads to the formation of a new residual vector until the tolerance is reached.

Once $\underline{\underline{\Theta}}_{n+1}$ has been obtained, the state laws allow the computation of the stress variables and the next step can be performed.

Finally, the model is implemented in the commercial finite element code Abaqus Explicit thanks to a user subroutine VUMAT dedicated to the material behaviour.

### 3 Parameters identification of the model

The identification procedure consists in determining, thanks to experimental data, the values of the mathematical constants involved in the constitutive relations of a model. To perform the identification, two major different methods depending on the complexity of the material and its associated model can be used. First, when the following constraining conditions can be met:
* possibility to associate a physical or phenomenological interpretation to a given parameter and a given test (Young modulus to account for linear elasticity during a uniaxial test...);

* possibility to make the distinction between the different physical phenomena (inelastic strain, loss of rigidity, rate dependence...);

* ability to make simplifying assumptions;

* capacity to use different experimental curves in a sequential order to find the parameters;

a sequential identification procedure can be developed. This method is quite empirical and gives very satisfying results for numerous models. However, in this context, the identification procedure is specific to a given material and its associated constitutive model and it is not possible to consider strongly coupled phenomena.

In our case, all the aspects of the behaviour are imbricated and it is not a simple matter to isolate the influence of each parameter in the global rate dependent response of the specimen. To circumvent this difficulty, one way is to consider the identification procedure as an optimization problem. It constitutes the second major class of methods.

The basic idea consists in minimizing the gap between experimental data and results coming from numerical simulations. In order to achieve this goal, the first step is to define the variables of comparison. For example, one can compare the experimental and numerical stress versus strain curves or, with recent optic measurement methods, fields of displacements. The next step is then to create a "cost function" that estimates the error made by the model while simulating the observed phenomena. Such an indicator can be an error within the meaning of least squares. Finally, from the knowledge of some sets of parameters and of the quality of their fit, a new set has to be determined in order to improve the numerical solution.

In this study, the identification procedure is based on a Pattern Search Algorithm which comes from the direct search (gradient-free) methods and enables to find minima of functions. This algorithm has been chosen for three principal qualities :

* its ability to deal with numerous variables;

* its robustness (if an incompatible sets of parameter

leading to non real cost result is encountered, the algorithm goes through the trap);

* and for its capacity to check a wide area in the parameters space which limits the starting point dependence in comparison with other methods such as descent ones.

Another point of interest in this algorithm is that the gradient of the cost function is not needed. However, it has to be noticed that the minimum found is a local one and so, the solution obtained is a particular solution that is expected to give a satisfying behaviour to the model.

In what follows, our attention is focused on the identification of the model in shear direction which is the part involving the most coefficients. Similar developments are used to identify the model in fibre direction for both tensile and compressive behaviours.

## 3.1 Experimental data available

In this first study "pseudo experimental data" are generated from the model described in [Rozycki (2000)]. This model includes both plasticity and static damage under the same framework as the present one. It traduces the mechanical behaviour of glass-fibre / epoxy-resin woven fabrics. A fundamental point in an identification procedure is to solicit all the aspects of the model as was shown by [Fossum (1997); Mahnken and Stein (1996)]. For composite woven fabrics, this involves irreversible strains, loss of stiffness, and rate dependence. To match these requirements, the generated "pseudo experimental data" represent cyclic shear tests for strain rates evolving

<figure><img src="image_2.png" /><figcaption>**Figure 1:** "Pseudo experimental data" for identification in shear direction</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Pattern Search Algorithm</figcaption></figure>

from $1 \cdot 10^{-3} s^{-1}$ to $1 \cdot 10^3 s^{-1}$. When the strain exceeds $0.1 \text{ m.m}^{-1}$, the strain rate is multiplied by a factor 10 in order to test more efficiently the strain rate effects. The stress versus strain curves obtained are plotted on figure (1).

The "pseudo experimental data" are then computed in tables with the time and the associated total strain and stress. For each given strain rate, a table named `data<sup>XP</sup>` = `[t; ε'; σ<sup>XP</sup>]<sub>ε</sub> is generated for comparison with numerical results.

## 3.2 Pattern Search Method

This section deals with the Pattern Search Algorithm. Let $\underline{\chi}$ be a set of parameters that is part of $\Upsilon$ the space of the admissible parameters delimited by lower and upper bounds. As presented below, $\Upsilon$ is composed of all the parameters involved in the shear direction of the model.

$$ \underline{\chi} = [G_{12}; v_{12}^G; \sigma^Y; K; m; Q; \beta; \alpha_{12}; \tau_{12}; Y_{12}^o; v_{12}^Y; Y_{12}^C; v_{12}^C; w_{12}^Y]^T $$

$$ \Upsilon = \{ \underline{\chi} / \chi_i^{\text{lower}} \le \chi_i \le \chi_i^{\text{upper}}, i \in \{1; \text{dim}(\underline{\chi})\} \} $$

For each set of parameters the time $t$ and the strain $\epsilon^t$ are used to drive the numerical simulations. Tables of numerical results are then computed under the form `dataNum` = `[t; ε'; σNum]`. It has to be noticed that the identification procedure is performed at the level of an integration point. By subtraction, stress error vectors $\underline{\Omega}$ are introduced and calculated by:

$$ \underline{\Omega} = \underline{\sigma}^{\text{XP}} - \underline{\sigma}^{\text{Num}} \qquad (18) $$

where each component of $\underline{\Omega}$ stands for the error on one point between the experimental and numerical values of stress. A global indicator of the quality of the identification is then computed from the norm of $\underline{\Omega}$. It takes into account all the experimental data available. In the least squares means, the indicator $||\underline{\Omega}||$ can be defined as:

$$ \Gamma = ||\underline{\Omega}|| = \frac{1}{2} \underline{\Omega} \cdot \underline{\Omega} \qquad (19) $$

To summarize, the problem can be expressed as: "Find $\underline{\chi}$ the set of parameters whose cost function $\Gamma$ is
<figure><img src="image_4.png" /><figcaption>Figure 3: Cost and associated number of iterations obtained for various initial sets of parameters</figcaption></figure>

minimum with $\underline{\chi} \in \Upsilon$ the space of the admissible parameters.”

To solve this problem, the Pattern Search Algorithm presented in figure (2) is employed.

Starting from a given set of parameters $\underline{\chi}^k$, the algorithm computes the cost and then prospects for a better candidate $\underline{\chi}^{k+1}$ by forming a "mesh" of points generated by the translation of $\underline{\chi}^k$ in each direction $e_i$ of the space of the parameters at a distance of $a_i^k$. If a better candidate appears, then a dilatation phase starts around $\underline{\chi}^{k+1}$ in the meaning that the next increment will be larger than the previous one ($||\underline{a}^{k+1}|| > ||\underline{a}^k||$). If there is no better candidate, the algorithm contracts its search by decreasing the distance of prospect ($||\underline{a}^{k+1}|| < ||\underline{a}^k||$). The procedure goes on until a tolerance on $||\underline{a}||$ or a maximum number of iterations is reached. For more information, consult [Lewis, Torczon, and Trosset (2000)].

### 3.3 Implementation of the method

The identification procedure is performed with the commercial software Matlab and the Genetic Algorithm Toolbox that includes a Pattern Search Algorithm. To conduct the procedure to success, some simplifications and assumptions have to be made on the set of parameters to identify. First of all, the huge vector of parameters:

$$ \underline{\chi} = [G_{12}; v_{12}^G; \sigma^y; K; m; Q; \beta; \alpha_{12}; \tau_{12}; Y_{12}^o; v_{12}^o; Y_{12}^c; v_{12}^c; w_{12}^c]^T $$

can be reduced with special considerations:

* $G_{12}$ the elastic modulus and $v_{12}^G$ its associated parameter that drives the strain rate dependence can be identified separately from the other parameters

by considering only the linear reversible part of the experimental data;

* $\tau_{12}$ that upper bounds the damage rate and $\alpha_{12}$ can be estimated "a priori" from the literature [Ladevèze, Allix, Deü, and Lévêque (2000), Suffis, Lubrecht, and Combescure (2003)] since there is, to our knowledge, no simple experimental test that enables to give a value to this parameters;

Finally, and with the preceding assumptions, the set of parameters to identify becomes:

$$ \underline{\chi} = [\sigma^y; K; m; Q; \beta; Y_{12}^o; v_{12}^o; Y_{12}^c; v_{12}^c; w_{12}^c]^T $$

A huge campaign of identification has been performed on several sets of data. The upper and lower bounds of the parameters were determined upon literature and physical considerations:

$$ \underline{\chi}^{lower} = [0; 0.1; 1.01; 0.1; 0.1; 0.001; 0.001; 0.1; 0.001; 0.001]^T $$

$$ \underline{\chi}^{upper} = [30; 2000; 5; 25000; 0.99; 1; 1; 1000; 1; 1]^T $$

The figure (3) shows the final set of parameters and its associated cost versus the position $\xi$ of the initial set of parameters $\underline{\chi}^0$ from the bounds. $\xi = 0$ means that $\underline{\chi}^0 = \underline{\chi}^{lower}$ and $\xi = 1$ means that $\underline{\chi}^0 = \underline{\chi}^{upper}$. The calculations are performed for $\xi = \{0; 0.25; 0.5; 0.75; 1\}$ with $G_{12} = 6396$ MPa, $v_{12}^G = 0.015$ (obtained from a first step of minimization upon the linear reversible parts of the experimental curves), $\tau = 1 \cdot 10^{-6} m^2 \cdot (m^{-2} \cdot s^{-1})$ and a reference strain rate $\dot{\epsilon}_{ref} = 1 \cdot 10^{-3} m \cdot (m \cdot s)^{-1}$.

First of all, it is visible that two curves differ from the others: $\xi = 0$ and $\xi = 0.75$. Next, for the initial sets of
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Comparison between experimental and numerical curves for both shear and fibre directions</figcaption></figure>

parameters located in the middle of the interval ($\xi = 0.25$ and $\xi = 0.5$), the final sets of parameters obtained after the optimization procedure are close to the ones from the others. Concerning the parameters of the delayed damage mesomodel, it appears that $Y_{12}^o$, $v_{12}^o$, $Y_{12}^c$, $v_{12}^c$ and $w_{12}^c$ all lie between their lower bound and 12 percent of their upper bound. The relative independence of these results from the starting point, and the satisfying global costs, lead us to think that the adequate parameters for the delayed damage mesomodel are close to the one determined in this first step of identification. The upper bounds of these parameters are then reduced to 15 percent of their initial value in the next steps of optimization.

Now, regarding the part of the model that deals with irreversible strains and strain rate dependence (viscoplasticity), the conclusions are more delicate to draw.

$\beta$ and $Q$, and if we do not consider the curves $\xi = 0.75$ and $\xi = 0$, the results are quite homogeneous and close to the identified values obtained in Rozycki (2000) that were used to generate the "pseudo experimental data". As $\beta$ and $Q$ drive the isotropic hardening phenomenon, the capability of the Pattern Search Algorithm to find the same values starting from very different points is quite interesting. For $K$, the Pattern Search Algorithm always leads nearby the lower bound. Except for $\xi = 0.75$, the same results are obtained for $\sigma^y$ which actually means that there is no elasticity, and that viscoplasticity occurs immediately... Last, $m$ seems to be limited by the imposed upper bound.

From the preceeding considerations, three conclusions can be drawn:
* the parameters of the delayed damage mesomodel are almost identified (even if the model is coupled, various sets for the "viscoplastic part" do not seem to affect them);

* the isotropic hardening parameters $\beta$ and $Q$ have been determined, and their value is similar to the one used to generate the "pseudo experimental curves";

* efforts have to be done to finalize the optimization process on the viscoplastic part of the model.

After this first step, the bounds of the almost identified parameters are severely reduced, and the Pattern Search Algorithm is launched again with the objective to specify the values of the parameters of the viscoplastic part of the model. Several intervals of admissible parameters are generated with higher upper bound for $m$ and lower upper bound for $K$.

## 3.4 Results

Finally, after a last step of optimization with a descent method to ensure that the minimum found is really a local minimum, the obtained results can be observed in the following figure (4) for shear direction (left side) and fibre direction (right side). The associated cost is equal to 34.

They clearly highlight the capacity of the optimization procedure to determine the parameters of the model. Moreover, the order of magnitude of the identified parameters seems to be physically reasonable (although the initial bounds were rather broad), which can give more confidence in the method. The Pattern Search Algorithm has demonstrated its ability to help for the identification of models dealing with strongly coupled physical phenomena. It has led to usable sets of parameters for the model in shear and fibre directions as summarized in tables (1) and (2).

The model presented in this paper allows to perform simulations on composite woven fabric structures upon six decades of strain rate.

## 4 Applications

In the two first parts of this section, the case of the dynamic extension of a bar is studied to check the ability of

<table><caption>**Table 1:** Parameters for the model in fibre direction</caption><thead><tr><th colspan="2">Elastic properties</th><th colspan="2">Damage properties</th></tr></thead><tbody><tr><td>E<sub>11</sub></td><td>30780MPa</td><td>Y<sup>o</sup><sub>11</sub></td><td>2.523 MPa</td></tr><tr><td>v<sup>E</sup><sub>11</sub></td><td>0.0211</td><td>v<sup>o</sup><sub>11</sub></td><td>0.041</td></tr><tr><td>ρ</td><td>1400 kg.m<sup>-3</sup></td><td>Y<sup>c</sup><sub>11</sub></td><td>275.38 MPa</td></tr><tr><td>ρ<sub>11</sub></td><td>3.64 · 10<sup>-4</sup></td><td>v<sup>l</sup><sub>11</sub></td><td>-0.018</td></tr><tr><td></td><td></td><td>v<sup>c</sup><sub>11</sub></td><td>2.07 · 10<sup>-4</sup></td></tr><tr><td></td><td></td><td>α<sub>11</sub></td><td>1</td></tr><tr><td></td><td></td><td>τ<sub>11</sub></td><td>1 · 10<sup>-5</sup> m<sup>2</sup>.(m<sup>-2</sup>.s<sup>-1</sup>)</td></tr></tbody></table>

<table><caption>**Table 2:** Parameters for the model in shear direction</caption><thead><tr><th colspan="2">Elastic properties</th><th colspan="2">Damage properties</th><th colspan="2">Viscoplastic properties</th></tr></thead><tbody><tr><td>G<sub>12</sub></td><td>6396 MPa</td><td>Y<sup>o</sup><sub>12</sub></td><td>8.8 · 10<sup>-5</sup> MPa</td><td>σ<sup>3</sup></td><td>16 MPa</td></tr><tr><td>v<sup>G</sup><sub>12</sub></td><td>1.5 · 10<sup>-2</sup></td><td>v<sup>o</sup><sub>12</sub></td><td>0.17</td><td>K</td><td>0.41</td></tr><tr><td>ρ</td><td>1400 kg.m<sup>-3</sup></td><td>Y<sup>c</sup><sub>12</sub></td><td>281 MPa</td><td>m</td><td>1.79</td></tr><tr><td></td><td></td><td>v<sup>l</sup><sub>12</sub></td><td>7.9 · 10<sup>-2</sup></td><td>Q</td><td>17315 MPa</td></tr><tr><td></td><td></td><td>v<sup>c</sup><sub>12</sub></td><td>6.47 · 10<sup>-3</sup></td><td>β</td><td>0.79</td></tr><tr><td></td><td></td><td>α<sub>12</sub></td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>τ<sub>12</sub></td><td>1 · 10<sup>-6</sup></td><td></td><td></td></tr></tbody></table>

the model to avoid strain localization and mesh dependence. This unidirectional example allows to differentiate the behaviour of the elementary ply in fibre direction (first part) and in shear direction (second part). In the two following parts, structural calculations are presented. The third part deals with the impact of carbon-fibre/epoxy-resin woven fabric plates, and the fourth one is related to the dynamic crushing of composite tubes.

### 4.1 Tensile simulations on a bar - behaviour in fibre direction

The bar considered in this study measures 100 mm length for a cross section area of 1 mm². One extremity is clamped whereas a constant velocity is imposed to the other one. Mesh dependence is investigated through the use of three different meshes involving 100, 200 and 400 linear truss elements. The equation of motion is integrated through the time with a standard explicit central difference scheme. The material behaviour is implemented in the finite element code Abaqus via a VUMAT subroutine. The material parameters employed for the model in fibre direction are given in table (1).

Many numerical simulations have been performed, and convergence of the material behaviour integration scheme was achieved in all cases with satisfying accuracy in few iterations. The upper graph of the figure (5) presents the profile of the damage obtained after
<figure><img src="image_6.png" /><figcaption>**Figure 5** : Damage along the bar for the model in fibre direction - influence of the mesh and of the strain rate</figcaption></figure>

2.5 · 10⁻⁵ s along the bar for different meshes and for a prescribed velocity of 5 · 10⁵ mm.s⁻¹ (which corresponds to an average macroscopic strain rate of 5000 m.(m.s)⁻¹). On the lower graph, the damage distribution along the bar for a mesh of 100 elements is plotted for various strain rates.

The first observation that can be made about the preceding curves is that the mesh objectivity of the results is ensured. Indeed, a modification of the mesh size of a factor 4 does not influence the response of the model in terms of damage distribution. Moreover, the curves obtained for various strain rates show that the fully damaged zone is not confined into a single element when the mesh is fine enough to capture the internal length associated with the delayed damage mesomodel and when the strain rate is high enough to make the damage propagate. It also

has to be added that the positivity of the dissipation was verified in all cases.

In conclusion, one can say that:

* the model is mesh objective;

* spurious localization phenomenon is avoided.

Hence, the introduction of strain rate dependent moduli in the formulation of the delayed damage mesomodel has not removed the interesting properties of the initial model but has, on the contrary, enlarged the domain of validity of the model in terms of strain rates.

It can be noted that the integration step time has to be chosen with care, which means that:

* the elastic wave should not propagate to more than
<figure><img src="image_7.png" /><figcaption>Figure 6: Damage along the bar for the model in shear direction - influence of the mesh and of the strain rate</figcaption></figure>

one element at a step time of integration $\Delta t < \frac{min(l_e)}{\sqrt{E_{max}/\rho}}$ (with $l_e$ the characteristic lengths of the elements of the model);

* the damage should not pass from 0 to 1 in one step time (which is traduced by the condition $\dot{l}_{max}\Delta t < 1$ or written differently: $\Delta t < \tau$).

Finally, these considerations lead to:

$$ \Delta t = min \left( \frac{min(l_e)}{\sqrt{\frac{E_{max}}{\rho}}} ; \tau \right) \quad (20) $$

and within these conditions, the Newton-Raphson procedure always converges in less than 5 or 6 iterations depending on the mesh size.

## 4.2 Tensile simulations on a bar - behaviour in shear direction

This part goes on with the preceding example to study the behaviour of the elementary ply submitted to shear loadings. The parameters used for the calculations are given in table (2). The distribution of the damage along the bar is plotted in figure (6). The influence of mesh size is presented for an imposed velocity taken equal to $5 \cdot 10^5 mm.s^{-1}$ on the upper part and to $5 \cdot 10^7 mm.s^{-1}$ on the lower part.

The coupling between viscoplasticity and delayed damage mesomodel disturbs neither the independence of the results from the mesh size nor the absence of localization phenomenon. Moreover, to check the influence of the step time on the response of the model, several calculations have been performed. The top curves of figure (7) below focus on the distribution of the damage along the bar for a given loading but for different step times.
<figure><img src="image_8.png" /><figcaption>**Figure 7** : Damage and cumulated plastic strain distribution along the bar - influence of the step time and of the mesh</figcaption></figure>

As a result, it can be seen that a modification of the step
time of a factor 100 does not modify the response of the
model. The integration of the constitutive laws of be-
haviour then appears to be sufficiently reliable to face
more sophisticated simulations.

Eventually, and in order to show the mesh objectivity of
the model on a "viscoplatic" variable, the graph in bot-
tom of the figure (7) plots the repartition of the cumulated
plastic strain among the bar.

In conclusion, the uniaxial example of the extension of
a bar has demonstrated the capacity of the model to
give similar results whatever the meshes and the step
times are. This fundamental intrinsic characteristic of the
model allows to predict the behaviour of structures inde-
pendently from the spatial and temporal discretizations.
To go further on the analysis of the model, two structural
applications are now presented.

4.3 Impact of a composite plate

Impact on layered composite plates involves many
damage phenomena [Smojver and Sori (2006)]. In this
study, the simulation of the impact of a composite plate
is inspired from the work of Johnson *et al*. Johnson,
Pickett, and Rozycki (2001) performed within the
German Aerospace Center (DLR).

A square plate, of thickness 4.6 mm and length 300 mm,
lying on a rigid square frame of length 250 mm, is
impacted normaly to its surface by a sphere of mass
21 kg launched at an initial velocity of 6280 mm.s⁻¹.
The layup sequence is [ 0 45 0 45 0 45 0 45 ] where each
ply is a woven fabric with the properties described in
tables (1) and (2). The Poisson's ratio v₁₂ is taken equal
to 0.13.

It is important to notice that the experiment was per-
formed on a carbon-fibre / epoxy-resin woven fabric
whereas the simulations are done with a glass-fibre /
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Damage state in the direction 1 of the fibres after impact and evolution of the reaction force during the impact</figcaption></figure>

epoxy-resin material. As a consequence, the model presented in this paper will be qualitatively checked to show its ability to reproduce observed physical phenomena.

Concerning the finite element model, the plate is meshed with the four nodes, six degrees of freedom linear multilayer shell elements. Both the impactor and the square frame are modelled as rigid solids. The contact between these two solids and the plate is considered to be a hard normal contact with no frictional sliding. All the degrees of freedom of the square frame are constrained whereas the impactor is only free to translate normally to the glass-fibre / epoxy-resin woven fabric plate. The calculations are performed up to a time of 20 ms so that the penetration of the impactor can be observed. The figures (8) below plot the state of damage on the upper and lower sides of the plate in fibre direction 1, and the evolution of reaction force through the time for two initial velocity of 6280 and 2330 mm.s⁻¹ (in which case, the impactor rebounds on the plate).

The state of damage presented here, for a time $t$ equal

to 3 ms, corresponds to the beginning of the failure of the plate as can be seen on the right curve with the abrupt decrease of the reaction force. The left figure clearly shows that the lower side of the plate is more damaged than the upper one. This numerical result is correlated with experimental observations. Moreover, as in the following figures (9), the global kinematics of damage is well predicted by the model. A macroscopic cross-shaped crack propagates under the impactor and finally causes the failure of the plate.

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Fully damaged zone after impact - experimental carbon-fibre / epoxy-resin on the left - numerical -glass-fibre / epoxy-resin on the right</figcaption></figure>

## 4.4 Axial dynamic crushing of a composite tube

This last structural example studies the dynamic crushing of a square composite thin-walled tube. Its interest lies in the preponderant sollicitation of the fibres in compres-
<figure><img src="image_11.png" /><figcaption>**Figure 10** : Simulation of the crushing of a composite tube - mid-length collapse mode</figcaption></figure>

sion. Qualitative comparisons are made with the exper-
imental results of [Mamalis, Manolakos, Ioannidis, and
Papapostolou (2005)].

The section of the tube measures 100 by 100 mm and
its height is also equal to 100 mm. The staking se-
quence is [0]14 with layers made of carbon-fibre/epoxy-
resin woven fabric. The total thickness is 3.5 mm. The
thin-walled tube lies on a rigid ground and is impacted
by a mass of 1 ton launched at an initial velocity of
5400 mm.s⁻¹. The modelling of this experiment is done
in the same manner as before. Several discretizations of

the tube have been checked, all leading to similar results.
Ground and impactor were modelled as rigid solids and
contact conditions were imposed on the tube. The fol-
lowing figure (10) presents the macrocrack due to the
failure of fibres in compression and the evolution of the
force applied by the tube on the impactor. The calcula-
tion was stopped just after failure.

The global collapse mode of the structure, a propagation
of a mid-length macroscopic crack, is in agreement with
the experimental observations. Moreover, the brittle fail-
ure of the thin-walled composite tube after a crushing
displacement of 2 mm is in a correct order of magnitude. These results, despite the fact they were obtained from a model identified on glass-fiber / epoxy-resin experimental data, are quite encouraging.

## 5 Conclusion

A model dedicated to glass-fibre / epoxy-resin woven fabrics was presented in this paper. The brittle failure of the fibres under tensile sollicitation, their nonlinear response due to the microbuckling in compression, the loss of stiffness and the irreversible strains observed in shear were all taken into account, and that, for a large range of strain rates. Moreover, special attention was paid to the construction of the model in order to avoid both the mesh dependency and the strain localization shortcomings that are often met in finite element simulations associated with softening material behaviours. After a successful step of identification based on a particular direct search method of optimization, calculations of uniaxial tensile tests on bars have shown the robustness of the model. Two structural applications have then been presented:

* the impact on a composite plate;

* and the crushing of a thin-walled tube.

In the first example, the model has predicted the kinematics of the failure of the plate with great accuracy. In the second example, the brittle rupture of the tube at mid-length was founded, in perfect agreement with the experimental observations.

## References

Chung, S. W.; Choi, Y. J.; Kim, S. J. (2003): Computational Simulation of Localized Damage by Finite Element Remeshing based on Bubble Packing Method. *CMES: Computer Modeling in Engineering & Sciences*, vol. 4, no. 6, pp. 707–718.

Fossum, A. F. (1997): Parameter estimation for an internal variable model using nonlinear optimization and analytical response sensitivities. *Journal of Engineering Materials and Technology*, vol. 119, pp. 337–345.

Gornet, L.; Marguet, S.; Marckmann, G. (2006): Numerical modelling of Nomex honeycomb cores: failure

and effective elastic properties. *European Conference on Computational Mechanics 2006*.

Herakovitch, C. T. (1998): *Mechanics of Fibrous Composites*. John Wiley and Sons, Inc.

Johnson, A.; Pickett, A.; Rozycki, P. (2001): Computational methods for predicting damage in composite structures. *Composites Science and Technology*, vol. 61, pp. 2183–2192.

Ladevèze, P.; Allix, O.; Deü, J.-F.; Lévêque, D. (2000): A mesomodel for localization and damage computation in laminates. *Computer Methods in Applied Mechanics and Engineering*, vol. 183, pp. 105–122.

Ladevèze, P.; Allix, O.; Gornet, L.; Lévêque, D.; Perret, L. (1998): A computation Damage Mechanics Approach for laminates : Identification and Comparison with Experimental Results. *Damage Mechanics in Engineering Materials*, pp. 481–502.

Ladevèze, P.; LeDantec, E. (1992): Damage modelling of the elementary ply for laminated composites. *Composites Science and Technology*, vol. 43, no. 3, pp. 257–267.

Lemaître, J.; Chaboche, J. (1990): *Mechanics of solid materials*. Cambridge University Press.

Lewis, R. M.; Torczon, V.; Trosset, M. W. (2000): Direct search methods: then and now. *Journal of Computational and Applied Mathematics*, vol. 124, pp. 191–207.

Mahnken, R.; Stein, E. (1996): Parameter identification for viscoplastic models based on analytical derivatives of a least-squares functional and stability investigations. *International Journal of Plasticity*, vol. 2, no. 4, pp. 451–479.

Mamalis, A.; Manolakos, D.; Ioannidis, M.; Papaportolou, D. (2005): On the response of thin-walled CFRP composite tubular components subjected to static and dynamic axial compressive loadings: experimental. *Composite Structures*, vol. 69, pp. 407–420.

Perzyna, P. (1966): Fundamental problems in viscoplasticity. *Advances in Applied Mechanics*, vol. 9, pp. 243–377.

Rozycki, P. (2000): Contribution au développement de lois de comportement pour matriaux composites soumis à l'impact. *PhD Thesis, LAMIH Valenciennes, France*.
**Simo, J.; Hugues, T. J. R.** (1998): *Computational Inelasticity*. Springer, New York.

**Smojver, I.; Sori, J.** (2006): *Numerical Modelling of Damage Response of Layered Composite Plates*. *CMC: Computers, Materials & Continua*, vol. 3, no. 1, pp. 13–24.

**Suffis, A.; Lubrecht, T. A. A.; Combescure, A.** (2003): *Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length*. *International Journal of Solids and Structures*, vol. 40, no. 13-14, pp. 3463–3476.