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

HAL
open science

# A Rate Dependent Constitutive Model for Carbon-Fiber Reinforced Plastic Woven Fabrics

Steven Marguet, Patrick Rozycki, Laurent Gornet

► To cite this version:

Steven Marguet, Patrick Rozycki, Laurent Gornet. A Rate Dependent Constitutive Model for Carbon-Fiber Reinforced Plastic Woven Fabrics. Mechanics of Advanced Materials and Structures, 2007, 14 (8), pp.619-631. 10.1080/15376490701672443. hal-01004935

HAL Id: hal-01004935

https://hal.science/hal-01004935

Submitted on 14 Dec 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Rate Dependent Constitutive Model for Carbon-Fiber Reinforced Plastic Woven Fabrics

S. Marguet, P. Rozycki, and L. Gornet

Ecole Centrale de Nantes, France

This paper deals with the modelling until rupture of composite structures made of carbon-fiber/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 behavior of the fibers under tensile solicitations, their nonlinear behavior in compression as well as the strongly nonlinear and irreversible behavior 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 procedure to this study, the models 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.

**Keywords** 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 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 behavior under dynamic loadings such as slamming or wave impacts then appears to be determinant for the design of safe structures.

Carbon-fiber/epoxy-resin woven fabrics composite materials have been studied for many years. Their linear orthotropic elastic behavior is now well known as can be seen in [1]. 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 [2]. 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 deterioration of the composite materials and eventually lead to the failure of the whole ship have to be dealt with.

Since the end of the 80s, important work has been 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 fibers in tension;

* the local buckling of the fibers 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 fibers 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 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 [3–5]. 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 the representation of 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 spatial discretization. To circumvent this issue, several strategies have been tested. Nonlocal models [6], gradient dependent models [7] or delayed damage mesomodels [4] are some of them. The aim is to introduce an internal length linked, for the delayed damage
mesomodel, with a maximal rate of damage. As a consequence, the damage cannot 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 [8]. 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 [9] to deal with the strain rate dependency of the glass-fiber/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 carbon-fiber/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 behavior 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 shown 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 (fibers and resin) and the structure. Moreover, and thanks to the thin aspect of most of the composite structures, the assumption of plane stress is made. 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 fiber 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 $\underline{\epsilon}^t$ is assumed:

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

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

To take into account the irreversible strains, the loss of rigidity, the different tensile and compressive behavior in fiber direction and the strain rate dependency, the thermodynamic potential is derived from the Helmholtz's volumic free energy $\psi$, and

written as follows:

$$ \begin{cases} 
2\psi (\underline{\epsilon}^e; \underline{d}; p; (\dot{\underline{\epsilon}}^e; t)) \\
= [\epsilon_{11}^e \cdot (1 - d_{11}) \cdot \Lambda_{11}^0 (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e] P^+ (\sigma_{11}^{t-\Delta t}) \\
+ [\epsilon_{11}^e \cdot (1 - d_{11}^c) \cdot \Lambda_{11}^0 (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e] P^- (\sigma_{11}^{t-\Delta t}) \\
+ \epsilon_{11}^e \cdot \Lambda_{12}^0 (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e \\
+ [\epsilon_{22}^e \cdot (1 - d_{22}) \cdot \Lambda_{22}^0 (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e] P^+ (\sigma_{22}^{t-\Delta t}) \\
+ [\epsilon_{22}^e \cdot (1 - d_{22}^c) \cdot \Lambda_{22}^0 (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e] P^- (\sigma_{22}^{t-\Delta t}) \\
+ \epsilon_{22}^e \cdot \Lambda_{21}^0 (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{11}^e \\
+ (2\epsilon_{12}^e) \cdot (1 - d_{12}) \cdot \Lambda_{1212}^0 (\dot{\epsilon}_{12}^e) \cdot (2\epsilon_{12}^e) + 2h(p)
\end{cases} \quad (2) $$

where:

* $\dot{\epsilon}^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 undamaged material which is supposed to be the initial one: $\Lambda_{11}^0(\dot{\epsilon}_{11}) = \frac{E_{11} \cdot V_{11}(\dot{\epsilon}_{11})}{1 - \nu_{12}\nu_{21}}$; $\Lambda_{22}^0(\dot{\epsilon}_{22}) = \frac{E_{22} \cdot V_{22}(\dot{\epsilon}_{22})}{1 - \nu_{12}\nu_{21}}$; $\Lambda_{12}^0(\dot{\epsilon}_{11}; \dot{\epsilon}_{22}) = \Lambda_{21}^0(\dot{\epsilon}_{11}; \dot{\epsilon}_{22}) = \frac{\nu_{12} \cdot E_{22} \cdot V_{12}(\dot{\epsilon}_{11}; \dot{\epsilon}_{22})}{1 - \nu_{12}\nu_{21}}$; $\Lambda_{1212}^0(\dot{\epsilon}_{1212}) = G_{12} \cdot V_{1212}(\dot{\epsilon}_{1212})$ with $E_{11}, E_{22}$ and $G_{12}$ the elastic moduli and $\nu_{12}, \nu_{21}$ the Poisson's moduli;

* $d_{11}$ and $d_{22}$ are damage variables that drive the brittle rupture of the fibers 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 fibers, at time $t$, they are explicitly expressed as a function of the stresses: $d_{ii}^c = \rho_{ii} \cdot \sigma_{ii}^{t-\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|+|\dot{\epsilon}_{22}^e|^{\nu_{11}+\nu_{22}}}{\dot{\epsilon}_{ref}^e}))$ 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 Eq. (2) leading to the expressions of $\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:

$$ \left\{ \begin{array}{l} \underline{\sigma} = \frac{\partial \psi}{\partial \epsilon^e} \\ \underline{\underline{Y}} = -\frac{\partial \psi}{\partial d} \\ R = \frac{\partial \psi}{\partial p} \end{array} \right. \qquad (3) $$

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

$$ \left\{ \begin{array}{l} \sigma_{11} = [(1 - d_{11}) \cdot \Lambda_{11}^0 (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e] \cdot P^+ (\sigma_{11}^{t-\delta t}) \\ \qquad + [(1 - d_{11}^e) \cdot \Lambda_{11}^0 (\dot{\epsilon}_{11}^e) \cdot \epsilon_{11}^e] \cdot P^- (\sigma_{11}^{t-\delta t}) \\ \qquad + \Lambda_{12}^0 (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e \\ \sigma_{22} = [(1 - d_{22}) \cdot \Lambda_{22}^0 (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e] \cdot P^+ (\sigma_{22}^{t-\delta t}) \\ \qquad + [(1 - d_{22}^e) \cdot \Lambda_{22}^0 (\dot{\epsilon}_{22}^e) \cdot \epsilon_{22}^e] \cdot P^- (\sigma_{22}^{t-\delta t}) \\ \qquad + \Lambda_{21}^0 (\dot{\epsilon}_{11}^e; \dot{\epsilon}_{22}^e) \cdot \epsilon_{11}^e \\ \sigma_{12} = (1 - d_{12}) \cdot \Lambda_{1212}^0 (\dot{\epsilon}_{12}^e) \cdot (2 \cdot \epsilon_{12}^e) \\ Y_{11} = \frac{1}{2} \cdot \Lambda_{11}^0 (\dot{\epsilon}_{11}^e) \cdot (\epsilon_{11}^e)^2 \\ Y_{22} = \frac{1}{2} \cdot \Lambda_{22}^0 (\dot{\epsilon}_{22}^e) \cdot (\epsilon_{22}^e)^2 \\ Y_{12} = 2 \cdot \Lambda_{1212}^0 (\dot{\epsilon}_{12}^e) \cdot (\epsilon_{12}^e - \epsilon_{12}^{vp})^2 \\ R = Q \cdot p^\beta \end{array} \right. \qquad (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 [10]. The following expression is adopted in order to deal with both the rate dependent irreversible effects and the degradations:

$$ \phi(\underline{\sigma}; \underline{Y}; R; (\underline{d})) = \phi^{vp}(\underline{\sigma}; \underline{Y}; R; (\underline{d})) + \phi^d(\underline{\sigma}; \underline{Y}; 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 fibers and resin.

Introducing the Perzyna's viscoplastic theory [11], 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^y \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 undamaged material to obtain the same state of deformation as the one obtained with $\sigma_{12}$ on the damaged material.

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

$$ \left\{ \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} \right. \qquad (8) $$

where:

$$ \dot{\lambda}^{vp} = \left\langle \frac{f}{K} \right\rangle_{+}^{m} \qquad (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 \le 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 stated:

$$ \left\{ \begin{array}{ll} \dot{d}_{ii} = \frac{1}{\tau_{ii}} \cdot (1 - e^{-\alpha_{ii} \cdot (z_{ii})_{+}}) & \text{if } d < d_{ii}^{c} ; 0 \text{ elsewhere,} \\ & \text{and } ii = \{11; 22; 12\} \end{array} \right. \qquad (10) $$

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 fiber direction, $z_{ii} = \frac{\sqrt{Y_{ii}^e - Y_{ii}^{c}}}{\sqrt{Y_{ii}^e - Y_{ii}^{c}}} - d_{ii}$ ($ii = 11, 22$) whereas in shear direction $z_{12} = \frac{\ln(Y_{eq}) - \ln(Y_{12}^e)}{\ln(Y_{12}^e) - \ln(Y_{12}^c)} - d_{12}$ with $Y_{ii}^e$ the yield thermodynamic force that initiates the damage and $Y_{ii}^c$ a critical thermodynamic force.
The evolution equations of the damage variables $\underline{d}$ define a damage mesomodel as was introduced by Ladevèze and his coworkers [4]. It is particularly interesting for its ability to avoid the strain localization phenomenon which results in mesh dependence of the results (see [8] 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}}^{vp} - \underline{\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{\epsilon}'_{n+1}$, 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 fiber directions;

2. integration of the coupled viscoplastic/delayed damage mesomodel in shear thanks to a modified Radial Return Mapping algorithm [12].

For the behavior in fiber 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}$ and $d'_{22}$ are updated and are used to check for the failure of the fibers in tension. Under compressive state, a simple criterion based on a critical stress drives the rupture.

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

$$\left\{ \begin{array}{l} 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} \end{array} \right. \quad (12)$$

where the evolution laws give the last right hand term of Eqs. (12). To condense these expressions, it is introduced the

vector of unknowns $\underline{\Theta} = (d; \epsilon^{vp}; p)^T$ which leads to:

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

As the expression of $\dot{\underline{\Theta}}$ is a non linear expression of $\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{\Theta}_{n+1}^k$, is supposed to be known (for the first iteration $k=0$, the approximation $\underline{\Theta}_{n+1}^0 = \underline{\Theta}_n$ is used). This hypothesis leads to the formation of a residual vector:

$$\underline{\underline{Res}}_{n+1}^k = \underline{\Theta}_{n+1}^k - \underline{\Theta}_n - \Delta t_{n+1} \dot{\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 Eq. (14) conducts to the following expression of the increment in strain variables:

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

with

$$\left[ \frac{\partial \underline{Res}}{\partial \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 \dot{\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 \dot{\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 \dot{\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{\Theta}_{n+1}^k = \underline{\Theta}_{n+1}^k + \delta \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{\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 behavior.

### 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 behavior 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. One way to circumvent this difficulty 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 mesurement 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 behavior 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 fiber direction for both tensile and compressive behaviors.

## 3.1. Experimental Data Available

In this first study and while waiting for the beginning of an imminent experimental campaign, “pseudo experimental data” are generated from the model described in [9]. This model includes both plasticity and static damage under the same framework as the present one. Here it is really important to notice that the “pseudo experimental data” follows the mechanical behavior of glass-fiber/epoxy-resin woven fabrics. So, in all that follows, the identification will also be made for this kind of materials. However, the method developed below is quite general and will be applied soon on new experimental data to deal with carbon-fiber/epoxy-resin woven fabrics. A fundamental point in an identification procedure is to solicit all the aspects of the model as was shown by [13, 14]. 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 from $1 \cdot 10^{-3} s^{-1}$ to $1 \cdot 10^3 s^{-1}$. When the strain exceeds $0.1 \ 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 $\overline{\underline{\underline{\mathit{data}}}}_{\epsilon}^{\mathit{XP}} = [t; \epsilon^t; \sigma^{\mathit{XP}}]_{\epsilon}$ is generated for comparison with numerical results.

## 3.2. Pattern Search Method

This section deals with the Pattern Search Algorithm. Let $\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

<figure><img src="image_3.png" /><figcaption>FIG. 1. “Pseudo experimental data” for identification in shear direction.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Cost and associated number of iterations obtained for various initial sets of parameters.</figcaption></figure>

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 “à priori” from the literature [4, 8] 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}^{\text{lower}} = [0; 0.1; 1.01; 0.1; 0.1; 0.001; 0.001; 0.1; 0.001; 0.001]^T \\ \underline{\chi}^{\text{upper}} = [30; 2000; 5; 25000; 0.99; 1; 1; 1000; 1; 1]^T $$

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}^{\text{lower}}$ and $\xi = 1$ means that $\underline{\chi}^0 = \underline{\chi}^{\text{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_{12} = 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 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. For $\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 [9] 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 preceding 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”;
<figure><img src="image_5.png" /><figcaption>FIG. 4. Comparison between experimental and numerical curves for both shear and fiber directions.</figcaption></figure>

* efforts have to be made 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 Figure 4 for shear direction (left side) and fiber 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 with the identification of models dealing with strongly coupled physical phenomena. It has led to usable sets of parameters for the model in shear and fiber directions as summarized in Tables 1 and 2 below.

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 the model to avoid strain localization and mesh dependence. This unidirectional example allows to differentiate the behavior of the elementary ply in fiber 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-fiber/epoxy-resin woven fabric plates, and the fourth is related to the dynamic crushing of composite tubes.

#### 4.1. Tensile Simulations on a Bar—Behavior in Fiber 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 on 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 behavior is implemented in the finite element code Abaqus via a VUMAT subroutine. The material parameters employed for the model in fiber direction are given in Table 1.

Many numerical simulations have been performed, and convergence of the material behavior integration scheme was achieved in all cases with satisfying accuracy in few iterations. The left side of the Figure 5 presents the profile of the damage obtained after 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 right side, the damage distribution along the bar for a mesh of 100 elements is plotted for various strain rates.

<table><caption>TABLE 1
Parameters for the model in fiber direction</caption><thead><tr><th></th><th>Elastic properties</th><th></th><th>Damage properties</th></tr></thead><tbody><tr><td>E<sub>11</sub></td><td>30780 MPa</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>Y<sup>o</sup></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>Q<sub>11</sub></td><td>3.64 · 10<sup>-4</sup></td><td>v<sup>Y<sup>c</sup></sup><sub>11</sub></td><td>-0.018</td></tr><tr><td></td><td></td><td>vv<sup>Y<sup>c</sup></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></th><th>Elastic properties</th><th>Damage properties</th><th>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> 8.8 &middot; 10<sup>-3</sup> MPa</td><td>&sigma;<sup>y</sup> 16 MPa</td></tr><tr><td>v<sup>G</sup><sub>12</sub></td><td>1.5 &middot; 10<sup>-2</sup></td><td>v<sup>o</sup><sub>12</sub> 0.17</td><td>K 0.41</td></tr><tr><td>&rho;</td><td>1400 kg.m<sup>-3</sup></td><td>Y<sup>c</sup><sub>12</sub> 281 MPa</td><td>m 1.79</td></tr><tr><td></td><td></td><td>v<sup>c</sup><sub>12</sub> 7.9 &middot; 10<sup>-2</sup></td><td>Q 17315 MPa</td></tr><tr><td></td><td></td><td>w<sup>c</sup><sub>12</sub> 6.47 &middot; 10<sup>-3</sup></td><td>&beta; 0.79</td></tr><tr><td></td><td></td><td>&alpha;<sub>12</sub> 1</td><td></td></tr><tr><td></td><td></td><td>&tau;<sub>12</sub> 1 &middot; 10<sup>-6</sup></td><td></td></tr></tbody></table>

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 its domain of validity 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 one element at a step time of integration $\Delta t < \frac{\min(l_e)}{\sqrt{\frac{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{d}_{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—Behavior in Shear Direction

This part goes on with the preceding example to study the behavior 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

<figure><img src="image_6.png" /><figcaption>FIG. 5. Damage along the bar for the model in fiber direction—influence of the mesh and of the strain rate.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Damage along the bar for the model in shear direction—influence of the mesh and of the strain rate.</figcaption></figure>

imposed velocity taken equal to $5 \cdot 10^5 \text{ mm.s}^{-1}$ on the left and to $5 \cdot 10^7 \text{ mm.s}^{-1}$ on the right.

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 left curves of Figure 7 below focus on the distribution of the damage along the bar for a given loading but for different step times. 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 behavior 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 “viscoplastic” variable, the right side of 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 behavior of structures independently from the spatial and temporal discretizations.

To go further on the analysis of the model, two structural applications are now presented.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Damage and cumulated plastic strain distribution along the bar—influence of the step time and of the mesh.</figcaption></figure>

### 4.3. Impact of a Composite Plate

The simulation of the impact of a composite plate is inspired from the work of Johnson et al. [16] performed within the German Aerospace Center (DLR). A square plate, of thickness $4.6 \text{ mm}$ and length $300 \text{ mm}$, lying on a rigid square frame of length $250 \text{ mm}$, is impacted normally to its surface by a sphere of mass $21 \text{ kg}$ launched at an initial velocity of $6280 \text{ mm.s}^{-1}$. 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 $\nu_{12}$ is taken equal to 0.13.

It is important to notice that the experiment was performed on a carbon-fiber/epoxy-resin woven fabric whereas the simulations are done with a glass-fiber/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.

In the finite element model, the plate is meshed with the multilayer and reduced integration shell elements S4R available in Abaqus Explicit. These elements are 4 nodes linear elements with 6 degrees of freedom per node. Several spatial discretizations, leading to similar results (for fine enough meshes), are tested. Both the impactor and the square frame are modelled as rigid solids. The contact beween 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
<figure><img src="image_9.png" /><figcaption>FIG. 8. Damage state in the direction 1 of the fibers after impact and evolution of the reaction force during the impact.</figcaption></figure>

whereas the impactor is only free to translate normally to the glass-fiber/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. Figures 8 below plot the state of damage on the upper and lower sides of the plate in fiber direction 1, and the evolution of reaction force through the time for two initial velocities 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 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.

### 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 fibers in compression. Qualitative comparisons are made with the experimental results of Mamalis et al. [17].

The section of the tube measures 100 by 100 mm and its height is also equal to 100 mm. The staking sequence is [0]₁₄ with layers made of carbon-fiber/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. Figure 10 presents the macrocrack due to the failure of fibers in compression on

<figure><img src="image_10.png" /><figcaption>FIG. 9. Fully damaged zone after impact—experimental carbon-fiber/epoxy-resin on the left—numerical-glass-fiber/epoxy-resin on the right.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Simulation of the crushing of a composite tube—mid-length collapse mode.</figcaption></figure>

the left, and the evolution of the force applied by the tube on the impactor on the right. The calculation 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 failure 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 carbon-fiber/epoxy-resin woven fabrics was presented in this paper. The brittle failure of the fibers under tensile loading, 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 behaviors. 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.

To recall with the design of huge sailing ships, having a model able to reproduce various degradation phenomena is essential to ensure *à priori* the integrity of the structure submitted to dynamic loadings. Moreover, as numerical models of these boats involve numerous degrees of freedom, it appears really interesting

to have reliable predictions whatever the meshes used for spatial discretization and the time steps of integration are. As encouraging results were obtained on the modelling of the glass-fiber/epoxy-resin woven fabrics, it can be said that the development of the model looks promising.

The next step of this study will be an experimental campaign on carbon-fiber/epoxy-resin woven fabrics. Other outlooks for this study deal with the modelling of the Nomex honeycomb core which connects the two composite skins in the panels used in nautical construction. First developments have been made in [2] concerning the linear elastic behavior until rupture of these structural materials. Further ones, dealing with the rate effects, will come next. A delaminating interface between the skins and the core will also be investigated in order to be more accurate on the design of oceanic ships subjected to dynamic loadings.

## REFERENCES

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

2. Gornet, L., Marckmann, G., and Ollier, G. "Interactions modèles expériences sur des âmes nids d'abeilles Nomex®: Application au design d'un voilier multicoque de course océanique," *Revue des composites et des matériaux avancés*, Editions Hermès (2006).

3. Ladevèze, P., and Le Dantec, E., "Damage modelling of the elementary ply for laminated composites," *Composites Science and Technology* 43(3), 257-267 (1992).

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

5. Ladevèze, P., Allix, O., Gornet, L., Lévêque, D., and Perret, L., "A computation damage mechanics approach for laminates: Identification and comparison with experimental results," *Damage Mechanics in Engineering Materials*, G. Voyiadjis, J. W. Wu, J. L. Chaboche (eds.), Elsevier, 481-502 (1998).

6. Bazant, Z. P., and Cabot, P. G., "Nonlocal continuum damage, localization instability and convergence," *ASME, International Journal of Applied Mechanics* 55, 287-293 (1988).

7. Zhang, H., Zhanga, X., and Chen, J.-S., "A new algorithm for numerical solution of dynamic elastic-plastic hardening and softening problems," *International Journal of Computers & Structures* 81(17), 1739-1749 (2003).
8. Suffis, A., Lubrecht, T. A. A., and Combescure, A., "Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length," *International Journal of Solids and Structures* **40**, 13–14, 3463–3476 (2003).

9. Rozycki, P., Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact, PhD Thesis, LAMIH Valenciennes, 2000.

10. Lemaître, J., and Chaboche, J. L., *Mechanics of solid materials*, Cambridge University Press (1990).

11. Perzyna, P., "Fundamental problems in viscoplasticity," *Advances in Applied Mechanics* **9**, 243–377 (1966).

12. Simo, J. C., and Hugues, T. J. R., *Computational inelasticity*, Springer, New York, 1998.

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

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

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

16. Johnson, A. F., Pickett, A. K., and Rozycki, P., "Computational methods for predicting damage in composite structures," *Composites Science and Technology* **61**, 2183–2192 (2001).

17. Mamalis, A. G., Manolakos, D. E., Ioannidis, M. B., and Papapostolou, D.P., "On the response of thin-walled CFRP composite tubular components subjected to static and dynamic axial compressive loadings: Experimental," *Composite Structures* **69**, 407–420 (2005).