# Different approaches for woven composite reinforcement forming simulation

Philippe Boisse, Nahiène Hamila, Fabrice Helenon, Benjamin Hagège, Jian Cao

► To cite this version:

Philippe Boisse, Nahiène Hamila, Fabrice Helenon, Benjamin Hagège, Jian Cao. Different approaches for woven composite reinforcement forming simulation. Ninth International ESAFORM conference on Material Forming, 2006, Glasgow, United Kingdom. pp.775-778, 10.1007/s12289-008-0002-7. hal-00081077

HAL Id: hal-00081077

https://hal.science/hal-00081077v1

Submitted on 15 Apr 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Different approaches for woven composite reinforcement forming simulation

P. Boisse • N. Hamila • F. Helenon • B. Hagege • J. Cao

**Abstract** Different approaches used for the simulation of woven reinforcement forming are investigated. Especially several methods based on finite element approximation are presented. Some are based on continuous modelling, while others, called discrete or mesoscopic approaches, model the components of the fabric. A semi discrete finite element made of woven unit cells under biaxial tension and in-plane shear is detailed. In continuous approaches, the difficulty lies in the necessity to take the strong specificity of the fibrous material into account. The yarn directions must be strictly followed during the large strains of the fabric. This is the main goal of the non-orthogonal model and of the hypoelastic constitutive model based on the yarn rotation presented in this paper. In the case of discrete and semi-discrete approaches the directions of the yarns are “naturally” followed because the yarns are modeled. Explicitly, however, modeling each component at the mesoscopic scale can lead to high numerical cost.

**Keywords** Fabrics/textiles · Finite element analysis · Composite forming · Continuous/discrete

P. Boisse (✉) · N. Hamila · F. Helenon
Laboratoire de Mécanique des Contacts et des Solides,
INSA de Lyon,
Villeurbanne Cedex, France
e-mail: Philippe.Boisse@insa-lyon.fr

J. Cao
Department of Mechanical Engineering, Northwestern University,
2145 Sheridan Road,
Evanston, IL 60208, USA

B. Hagege
Laboratoire de Mécanique Roberval,
Université Technologique de Compiègne,
Compiegne, France

## Introduction

In order to determine the deformed shape of draped fabrics, several codes have been developed based on geometrical approaches so called fishnet algorithms [1–4]. These methods, where the fabric is placed progressively from an initial line, provide a close enough resemblance to handmade draping. They are very fast and fairly efficient in many prepreg draping cases. Nevertheless, these methods have major drawbacks. They account neither for the mechanical behaviour of the fabric nor for the static boundary conditions. This last point is very important in the case of forming with punch and die (such as in the preforming of the R.T.M. process). The loads on the tools, especially on the blank-holder, influence the quality of the shaping operation, and therefore, need to be considered in simulations.

The alternative to these geometrical methods consists of a mechanical analysis of the fabric deformation under the boundary conditions prescribed by the forming process. This requires a model of the woven reinforcement and its mechanical behaviour, in order to achieve the deformation through a numerical method, for instance, the Finite Element Method. The mechanical behaviour of fabrics is complex due to the intricate interactions of the yarns. It is a multi-scale problem. The macroscopic behaviour is very much dependent on the interactions of yarns at the meso-scale (scale of the woven unit cell) and at the micro-scale (level of the fibres constituting yarns). Despite of a great amount of work in the field, there is no widely accepted model that accurately describes all the main aspects of fabric mechanical behaviour. The main model families come from the multi-scale nature of the textile. A first family of models is obtained by homogenizing the mechanical behaviour of the underlying meso-structure and considering the fabric as an anisotropic continuum
[5–13]. If these models can easily be integrated in standard finite element (FE) shell or membrane elements, then the identification of homogenized material parameters is difficult, especially because these parameters change when the fabric is strained and when, consequently, the directions and the geometry (crimp, transverse sections...) of the yarns change. Some of these approaches will be described, especially a non-orthogonal constitutive model [8, 9] and an anisotropic hypoelastic continuous behaviour for fibrous material based on an objective derivative using the rotation of the fibre [10, 11].

Conversely, some authors present fully discrete models of fabrics [14–18]. Each yarn or each fibre is modelled and is assumed to be a straight or a curved beam or truss. Sometimes they are modeled as 3D domains [16]. Springs are often used to model warp and weft yarn interactions. In the objective of fabric forming simulations, some authors extend the discrete modelling to the whole textile structure that is represented by a network of interwoven trusses or beams with different tensional and rotational springs. Accounting for the simplicity of each component, the whole textile structure deformation can be computed. Nevertheless, the computational effort needed is relatively significant. At present this method is restricted to simple geometry of the local yarn and relatively simple mechanical behaviour. When a fine model of the fibrous yarns is used, the analysis can only consider a small part of the textile reinforcement such as a few woven or knitted cells.

The semi-discrete approach is a compromise between the above continuous and discrete approaches [19, 20]. A finite element method is associated to a mesoscopic analysis of the woven unit cell. Specific finite elements are defined that are made of a discrete number of woven unit cells. The mechanical behaviour of these woven cells is obtained by experimental analyses or from 3D FE computations of the woven cell. The nodal interior loads are deduced from this local behaviour and the corresponding strain energy in the element deformation. One objective of the present paper is, to assess the importance of in-plane shear strain energy in fabric forming simulations. It will be shown that the in-plane shear stiffness is particularly important for the description of wrinkles and shear strain after the shear locking angle.

## Continuous approach for composite forming process analysis

F.E. analysis of composite forming requires modelling of the different aspects involved in the process and especially a constitutive mechanical model of the fibrous reinforcement. The multiscale nature of the composite and of its fibrous reinforcement permits different possible

approaches. The first one considers the fibrous reinforcement as a continuum. The reinforcement is not continuous at lower scales but it is usually continuous in average and a continuous material superimposed to the fibrous material can be considered. The advantage of the continuous approach is that it can be used in standard finite elements. Nevertheless, the constitutive model of this continuum will have to convey the very specific mechanical behaviour of the fibrous reinforcement. Especially this behaviour is mainly depending on the fibre directions that are strongly changing during forming. It has been shown that the approaches usually used for anisotropic metal forming cannot be used here [10]. Other specific aspects of fibrous reinforcement mechanical, behaviour such as crimp change and shear locking should also be taken into account in the model.

The main difficulty in using the continuous approach is capturing the effects of the fibre architecture and its evolution. There are many models. Most of them assume that the fibrous reinforcement is elastic while forming. That is usually true for extensions in the fibre directions, but not obvious in the other directions, such as in-plane shear, bending, and transverse compression. Nevertheless, the forming process is a more or less monotonous operation and making this assumption doesn't change the result of the analysis greatly. Continuous behaviour models that capture macro-level phenomena at lower scales generally concern homogenisation, although homogenisation typically refers to techniques applied to a two-scale periodic material, in which the analysis of a unit cell reveals the properties of the homogenised material [5]. This approach is elegant but requires lengthy computational times. Furthermore, extending it to non-linear problems is difficult. Two continuous approaches used in FE analysis are described below.

## Non-orthogonal constitutive models

In these approaches, the stress and strain of a continuous material are related to fibrous reinforcement using the constitutive relation in a non-orthogonal frame directed by the fibre directions. Models have been developed by [7–9]. They considered two yarn directions and used them to define the non-orthogonal frame. The model developed in [9] is briefly described here.

This approach uses the Green Naghdi frame. It is an orthonormal frame which is rotated by $R$, the rotation of the polar decomposition in which the local stress increment computations at finite strain are made. It is used by many F. E. codes especially ABAQUS. This basis is noted ($e_1, e_2$) (Fig. 1). A second non-orthogonal basis ($g_1, g_2$) is defined from the yarn directions. The contravariant basis ($g^1, g^2$) is associated ($g_\beta \cdot g^\alpha = \delta^\alpha_\beta$ ($\alpha, \beta=1$ or 2)). The ($e_1, e_2$) basis
<figure><img src="image_2.png" /><figcaption>Fig. 1 Schematic of undeformed and deformed fabric [9]</figcaption></figure>

permits an objective sum of the stress components in this frame:

$$[\sigma^{t+\bar{d}'}] = [\sigma^t] + [d\sigma] \quad (1)$$

The stress sum (Eq. 1) is made in the ($\mathbf{e}_1$, $\mathbf{e}_2$) basis, but the mechanical behaviour of the woven reinforcement is

$$[T_1] = \begin{bmatrix} \cos^2 \alpha & \sin^2 \alpha & \sin \alpha \cos \alpha \\ \cos^2 (\alpha + \theta) & \sin^2 (\alpha + \theta) & \sin (\alpha + \theta) \cos (\alpha + \theta) \\ 2 \cos \alpha \cos (\alpha + \theta) & 2 \sin \alpha \sin (\alpha + \theta) & \sin (2\alpha + \theta) \end{bmatrix} \quad (5)$$

Consequently:

$$[d\sigma] = [T_1]^T [\tilde{D}] [T_1] [d\epsilon] = [D] [d\epsilon] \quad (6)$$

$[D]$ is the constitutive matrix in ($\mathbf{e}_1$, $\mathbf{e}_2$) basis and the new value of the stress can be computed by Eq. 1. The calculation of the $[T_1]$ matrix needs the knowledge of $\alpha$ and $\theta$ angles (Fig. 1). Those are calculated respectively from $\mathbf{F}$ (the deformation gradient tensor) and $\mathbf{R}$ (the polar rotation):

$$\mathbf{e}_i^t = \mathbf{R} \cdot \mathbf{e}_i^0 \quad \mathbf{g}_i^t = \mathbf{F} \cdot \mathbf{g}_i^0 \quad (7)$$

Based on this approach, the simulation of a picture frame shearing of a woven composite fabric is shown Fig. 2b and c [9]. It is compared to experimental results [21] (Fig. 2).

An hypoelastic model for fibrous materials

Approaches traditionally developed in finite element codes for anisotropic metal at large strains are based on Jaumann corotational formulation [22–24] or the Green–Naghdi approach [24–26]. In these models, a rotation is used both to define an objective derivative for the hypoelastic law and to update the orthotropic frame. The rotations used in Green–Naghdi and Jaumann derivatives are average rotations of the material (polar rotation and corotational rotations respectively). These well-known approaches cannot be used for a fibrous material under large strains because the update for the strong direction must follow the

known in the frame of the yarns ($\mathbf{g}_1$, $\mathbf{g}_2$). The following sets of components are considered:

$$\begin{align}
\sigma &= \sigma_{\alpha\beta} \mathbf{e}_\alpha \otimes \mathbf{e}_\beta = \tilde{\sigma}^{\alpha\beta} \mathbf{g}_\alpha \otimes \mathbf{g}_\beta \nonumber \\
\epsilon &= \epsilon_{\alpha\beta} \mathbf{e}_\alpha \otimes \mathbf{e}_\beta = \tilde{\epsilon}^{\alpha\beta} \mathbf{g}_\alpha \otimes \mathbf{g}_\beta \tag{2}
\end{align}$$

Assuming an orthotropic format of the constitutive matrix in the non-orthogonal frame,

$$\begin{bmatrix} d\tilde{\sigma}^{11} \\ d\tilde{\sigma}^{22} \\ d\tilde{\sigma}^{12} \end{bmatrix} = \begin{bmatrix} \tilde{D}^{11} & \tilde{D}^{12} & 0 \\ \tilde{D}^{12} & \tilde{D}^{22} & 0 \\ 0 & 0 & \tilde{D}^{33} \end{bmatrix} \begin{bmatrix} d\tilde{\epsilon}^{11} \\ d\tilde{\epsilon}^{22} \\ d\tilde{\epsilon}^{12} \end{bmatrix} \quad (3)$$

This relation can be written in the ($\mathbf{e}_1$, $\mathbf{e}_2$) basis. Equation 2 leads to:

$$[d\tilde{\epsilon}] = [T_1][d\epsilon] \quad [d\sigma] = [T_1]^T [d\tilde{\sigma}] \quad (4)$$

where $[T_1]$ is the following matrix (given by Eq. 2):

fibre direction strictly [10]. The following approach uses the rotation of the fibre denoted $\Delta$ [10, 11]. The rotation $\Delta$ is used to update the initial constitutive axes {$\mathbf{k}^0$} to the current constitutive axes {$\mathbf{k}^t$}

$$\mathbf{k}_i^t = \Delta \cdot \mathbf{k}_i^0 \quad (8)$$

Equation 8 leads to equations (Eq. 9) that explicitly give the constitutive axes {$\mathbf{k}^t$} as functions of the initial constitutive axes {$\mathbf{k}^0$} and the deformation gradient $\mathbf{F}$ [24]:

$$\begin{align}
\mathbf{k}_1^t &= \frac{\mathbf{F} \cdot \mathbf{k}_1^0}{\|\mathbf{F} \cdot \mathbf{k}_1^0\|} \qquad \mathbf{k}_2^t = \mathbf{k}_2^0 - \frac{b_2}{1+b_1} (\mathbf{k}_1^0 + \mathbf{k}_1^t) \nonumber \\
\mathbf{k}_3^t &= \mathbf{k}_3^0 - \frac{b_3}{1+b_1} (\mathbf{k}_1^0 + \mathbf{k}_1^t) \tag{9}
\end{align}$$

with $b_k = \mathbf{k}_1^t \cdot \mathbf{k}_k^0$ and $b_k \neq 1$. In this formulation, the fibre direction, i.e. the strong anisotropic direction, remains aligned with the first vector of {$\mathbf{k}^t$}. The constitutive behaviour is then fully defined at each time point. In fact, the component of the initial constitutive tensor $^0\mathbf{C}$ can be computed from the traditional engineer's constants:

$$^0\mathbf{C} = {}^0C_{ijkl} \mathbf{k}_i^0 \otimes \mathbf{k}_j^0 \otimes \mathbf{k}_k^0 \otimes \mathbf{k}_l^0 \quad (10)$$

The current constitutive tensor $\mathbf{C}$ can be deduced from $^0\mathbf{C}$ by a rotational transport based on the fourth order rotation tensor $\Delta$:

$$\mathbf{C} = \Lambda : {}^0\mathbf{C} : \Lambda^T \quad \text{with} \quad \forall \mathbf{A} \ \Lambda : \mathbf{A} = \Delta \cdot \mathbf{A} \cdot \Delta^T \quad (11)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2 a Picture frame in initial and deformed state [21]. b FEM model for trellising composite fabrics under shear frame. c Test/FEM comparison [9]</figcaption></figure>

Consequently C is given by:

$$C = {}^{0}C_{ijkl}k_{i}^{t} \otimes k_{j}^{t} \otimes k_{k}^{t} \otimes k_{l}^{t} \qquad (12)$$

The constitutive tensor C is used in a hypo-elastic law written:

$$\sigma^{\nabla} = C : D \quad \text{with} \quad \sigma^{\nabla} = \Delta \cdot \frac{d}{dt}(\Delta^{T} \cdot \sigma \cdot \Delta) \cdot \Delta^{T} \qquad (13)$$

D is the strain rate and $\sigma^{\nabla}$ is the objective derivative of the Cauchy stress associated to the fibre rotation $\Delta$. The

<figure><img src="image_4.png" /><figcaption>Fig. 3 Simulation result of a bias extension test on a NCF using a mesoscopic or discrete approach [16]</figcaption></figure>

cumulated tensorial strain tensor $\varepsilon$ and stress tensor $\sigma$ associated with such an objective derivative are:

$$\varepsilon = \Delta \cdot (\int_{0}^{t} \Delta^{T} \cdot D \cdot \Delta dt) \cdot \Delta^{T} \qquad (14)$$

$$\sigma = \Delta \cdot (\int_{0}^{t} \Delta^{T} \cdot C \cdot D \cdot \Delta dt) \cdot \Delta^{T}$$

It can be shown that Eq. 14 will always give a logarithmic strain in the strong anisotropic direction and that it ensures the summation of the stress increments along this direction [27]. Finally, the use of the material rotation tensor $\Delta$ for the
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Four node finite element made of woven cells</figcaption></figure>

objective derivative (Eq. 13) and the evolution law (Eq. 11) entails a consistent approach for fibrous media.

The above approach is carried out for a single fibre direction, although many fibrous reinforcements, and especially woven fabrics, have two fibre directions. The behaviour here is no longer orthotropic because there are very large angle variations between the warp and weft directions due to in-plane shear. The above formulation can be used by superposition on the same point for two materials using their own fibre directions. Sliding between warp and weft directions can be ignored because the two materials (warp and weft) are in the same finite element. However, that does not take account of the interaction between fibres, such as crimp changes.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Deep drawing of a square box: geometry of the tools</figcaption></figure>

## Discrete or mesoscopic approach

The discrete approach is the opposite of the continuous approach, it considers and models the components of fibrous reinforcement at a low scale. These components can be yarns, woven cells or stitching, and also sometimes fibres. Because these elements are usually at the mesoscale, the approach is also known as meso-mechanical modelling. Some analyses have been proposed where all fibres are modelled (microscale modelling) [14, 17] but the number of fibres in a composite structure limits these computations to small sub-domains, for instance a woven cell or a few braided or knitted loops. Discrete approaches for textile forming are generally made at mesoscopic scale [15, 16, 18]. A major difficulty lies in the description of the components at mesoscopic scale, usually the woven yarns. A compromise must be found between a precise description (which will be expensive from the computation time point of view) and a simple description, where it is possible to compute the entire forming process.

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Deformed shape and shear angle (tensile energy only)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Deformed shape and shear angle (tensile + in plane shear energy)</figcaption></figure>

Beam, truss or spring elements are the more common descriptions for yarns [14, 15, 17, 18]. In [16], mesomechanical FE modelling of NCF (Non Crimp fabric) uses 3D elements for each yarn and bar element for stitching (Fig. 3). Friction contact between tows and plies is taken into account. The simulation of a bias test is shown in Fig. 3. The simulation of this test has proved to be difficult in case of continuous approaches [28, 29]. Although the complete model probably needs extensive computational time, a forming process has been simulated in this way. A point in favour of analysis at the mesoscopic scale lies in the strong increase in computer efficiency.

## Semi discrete finite element

The approach associates a finite element method with a mesoscopic analysis of the woven unit cell. Specific finite elements are defined that are made of a discrete number of woven unit cells. The description of the fabric by finite elements needs to assume that two points of a weft and a warp yarns initially superimposed remain superimposed after forming, i.e., there is no translation sliding between the yarns. That has been experimentally shown in most cases [30]. Nevertheless, the model will have to describe the specificities of textile reinforcement mechanical behaviour, especially:

* the non-linear tensile behaviour due to crimp interchange

* shear locking angle and the very different in plane shear behaviour before and after this angle.

A global simplified dynamic equation is considered:

$$ 
\begin{aligned}
& \sum_{p=1}^{\text{ncell}} (p\varepsilon_{11}(\boldsymbol{\eta})^P T^{11} L_1 + p\varepsilon_{22}(\boldsymbol{\eta})^P T^{22} L_2) \\
& + \sum_{p=1}^{\text{ncell}} p C^p \gamma(\boldsymbol{\eta}) - W_{\text{ext}}(\boldsymbol{\eta}) = -\int_{\Omega} \rho \ddot{\boldsymbol{u}} \cdot \boldsymbol{\eta} \, dV \quad (15) \\
\forall \boldsymbol{\eta} \text{ such as } \boldsymbol{\eta}=0 \text{ on } \Gamma_u
\end{aligned}
$$

where $\varepsilon(\boldsymbol{\eta}) = \nabla^s \boldsymbol{\eta} = \varepsilon_{\alpha\beta}(\boldsymbol{\eta})\boldsymbol{h}^\alpha \otimes \boldsymbol{h}^\beta$ is the symmetrical gradient in the virtual displacement $\boldsymbol{\eta}$ ($\alpha$ and $\beta$ are indexes taking value 1 or 2). $\boldsymbol{h}^1$, $\boldsymbol{h}^2$ are the contravariant vectors related to the unit vector in warp and weft directions $\boldsymbol{h}_1$ and $\boldsymbol{h}_2$. $L_1$ and $L_2$ are the lengths of the warp and weft yarns in the mid-plane of the fabric. $\gamma(\boldsymbol{\eta})$ is the virtual relative rotation between warp and weft fibres (or virtual shear angle). ncell is the number of woven unit cells of the textile structure, $^pQ$, means that the quantity $Q$ is considered for the woven unit cell number $p$. $\rho$ is the mass per volume of the fabric $\Omega$. $W_{\text{ext}}(\boldsymbol{\eta})$ is the virtual work of the exterior prescribed loads. $\Gamma_u$ is the boundary with prescribed displacements.

To make a finite element simulation of composite woven reinforcement forming based on Eq. 15, it is necessary to be able to know the tensions $T^{11}$ and $T^{22}$ and the shear couple $C$ for a given strain field in the woven unit cell. It is assumed that the tension do not depend on the shear angle and that the shear couple do not depend on the axial strain i.e. $T^{11}$ ($\varepsilon_{11}$, $\varepsilon_{22}$), $T^{22}$ ($\varepsilon_{11}$, $\varepsilon_{22}$) and $C(\gamma)$. In [31], biaxial tensile tests performed for different angle between warp and weft yarns have shown that the influence of this angle is small and can be neglected. The second assumption ($C$ only depending on $\gamma$) is probably less true [32, 33]. Nevertheless, all the experimental results that are currently available give the shear load in function of the shear angle without any information on the tensions, so the assumption $C(\gamma)$ is made by default.

The tension surfaces ($T^{11}$ ($\varepsilon_{11}$, $\varepsilon_{22}$) and $T^{22}$ ($\varepsilon_{11}$, $\varepsilon_{22}$)) and of the shear curve ($C(\gamma)$) can be determined both by experiments or mesoscopic F.E. on the unit woven cell [11, 34]. From the simplified dynamic Eq. 15 specific finite elements for fabric forming are constructed (Fig. 4 in the case of a four-node quadrangle). These elements are composed of woven cells the tensile and in-plane shear strains of which are given by nodal displacements. From

<table><caption>**Table 1** Mechanical properties of the glass plain weave ("Deep drawing of a square box. Influence of the in plane shear stiffness" section)</caption><tr><td>Tensile stiffness in warp direction the tension is assumed to be</td><td>T<sup>11</sup>=C<sub>1</sub>ε<sub>11</sub></td><td>C<sub>1</sub> is assumed to be constant: C<sub>1</sub>=46,000 N/yarn</td></tr><tr><td>Tensile stiffness in weft direction the tension is assumed to be</td><td>T<sup>22</sup>=C<sub>2</sub>ε<sub>22</sub></td><td>C<sub>2</sub> is assumed to be constant: C<sub>2</sub>=48,000 N/yarn</td></tr><tr><td>In plane shear stiffness:</td><td>C(γ) = k<sub>1</sub> · γ si γ ≤ γ<sub>c</sub><br/>C(γ) = k<sub>2</sub> · (γ - γ<sub>c</sub>) + k<sub>1</sub> · γ<sub>c</sub> if γ > γ<sub>c</sub></td><td>k<sub>1</sub>=0.25 mm N/rad, k<sub>2</sub>=4.08 mm N/rad and, γ<sub>c</sub>=0.72 rad (41°)</td></tr></table>
**Fig. 8** Hemispherical forming of an unbalanced composite fabric. **a** Geometry of the tools.

**b** Final shape obtained by a simulation with semi-discrete finite elements. **c** Experimental final shape

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

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

Eq. 5 the nodal interior load components ($F_{\text{int}}^e$)$_s$ are obtained from the strain interpolation matrix:

$$W_{\text{int}}^{e}(\boldsymbol{\eta}) = \eta_{s} \sum_{p=1}^{\text{ncell}^{e}} \left( \|\boldsymbol{g}_{1}\|^{-2} p B_{11s} p T^{11} + \|\boldsymbol{g}_{2}\|^{-2} p B_{22s} p T^{22} + p C^{p} B_{\gamma s} \right) = \eta_{s} (F_{\text{int}}^{e})_{s} \quad (16)$$

These nodal interior loads are directly calculated in an explicit dynamic approach. The detailed formulation of the element can be found in [20]. An extension to a three-node element where the warp and weft directions of the woven fabric can be in arbitrary direction with respect to the direction of the element side is presented in [35].

Deep drawing of a square box. Influence of the in plane shear stiffness

This classical benchmark for sheet metal forming [36] (Fig. 5) is performed using the semi-discrete finite elements presented in “Semi discrete finite element” section. The geometry (Fig. 4 for a quarter of the system for symmetry reasons) is strongly double curved, and this test is severe especially for fabric forming because it results in very large angle variations between warp and weft yarns in the radius of the square box. The shear angles that are necessary to shape the part can be larger than the locking angle of the fabric depending on the radius values. The forming of a glass plain weave fabric is simulated. The mechanical

<table><caption>**Table 2** Mechanical properties of the unbalanced fabric (“Hemispherical forming of an unbalanced composite fabric” section)</caption><tr><td>Tensile stiffness in warp direction the tension is assumed to be</td><td>T<sup>11</sup>=C<sub>1</sub>ε<sub>11</sub></td><td>C<sub>1</sub> is assumed to be constant: C<sub>1</sub>=50 N/yarn</td></tr><tr><td>Tensile stiffness in weft direction the tension is assumed to be</td><td>T<sup>22</sup>=C<sub>2</sub>ε<sub>22</sub></td><td>C<sub>2</sub> is assumed to be constant: C<sub>2</sub>=0.2 N/yarn</td></tr><tr><td>In plane shear stiffness:</td><td>C(γ) = k<sub>1</sub>γ si γ ≤ γ<sub>c</sub><br>C(γ) = k<sub>2</sub> (γ - γ<sub>c</sub>) + k<sub>1</sub>γ<sub>c</sub> if γ > γ<sub>c</sub></td><td>k<sub>1</sub>=0.03 mm N/rad, k<sub>2</sub>=0.095 mm N/rad and, γ<sub>c</sub>=0.5 rad (29°)</td></tr></table>
properties are given in Table 1 [37]. The locking angle of
this fabric is 41°. Beyond this angle, the in-plane shear
stiffness becomes large. The forming simulation is made
with the dynamic equation presented in Eq. 15 (tensile and
in plane shear energy) but also neglecting the in-plane shear
term (tensile energy only). Figures 6 and 7 present in both
cases the computed deformed shape after forming. Because
of the geometry (strongly double curved) the locking angle
of the plain weave fabric is exceeded. This leads to rather
different results for both approaches. There is no wrinkle
when using the approach base on the only tension because
there is no source of instability (Fig. 6). On the contrary, the
computed solution obtained when taking shear into account
shows some wrinkles (Fig. 7). Those are due to shear
locking that leads to out-of plane solutions in order to
reduce this shear. The rotation angles are clearly reduced
when the shear stiffness is taken into account. This example
(as others [19]) shows that the contribution of the in-plane
shear term is principally in the description of the deformed
shape after the appearance of wrinkles (Fig. 7).

Hemispherical forming of an unbalanced composite fabric

The simulation of a very unbalanced fabric is presented
Fig. 8 (The simulation is performed on a quarter of the part
for symmetry reasons.). The warp yarns are much stiffer
than the weft yarns (Table 2 [37]). The blank holder is a
6 kg ring submitted to its own weight. The experimental
forming has been performed at Nottingham University [37].
The final shape obtained after forming is very asymmetri-
cal. There is a large axial strain in the weft direction
(horizontal) and large displacements with very small axial
strain in the warp direction (vertical). This final shape is
well obtained by the simulation. The ratio of the lengths
after deformation l<sub>wen</sub>/l<sub>warp</sub> is equal at the top of the
hemisphere to 1.8 in experiments and in simulation as well.
There are many wrinkles, especially along the vertical axis.
They are fairly well obtained by the simulation.

**Conclusion**

The simulation of composite reinforcement forming is a field
which undergoes many developments. Several approaches
exist and mainly differ according to the scale at which the
modelling is made. In continuous approaches, it is assumed
that a continuous medium can stand in for the set of yarns or
fibres. In this case, the difficulty lies in the necessity for the
macro model to take the strong specificity of the fibrous
material into account, especially that the yarn directions must
be strictly followed during the large strains of the fabric. This
is the main goal of the non-orthogonal model and of the
hypoelastic constitutive model based on the yarn rotation

presented in this paper. In the case of discrete and semi-
discrete approaches the directions of the yarns are “natural-
ly” followed because the yarns are themselves modelled. In
each approach different assumption and different level of
modelling can be done. It has been shown, for instance, that
the in-plane shear rigidity mainly influences the wrinkles
that appear when the shear angle exceeds the locking angle.

**Acknowledgements** The work reported here has been carried out in the scope of the project ITOOL (European Commission) and the National Science Foundation of US.

**References**

1. Mark C, Taylor HM (1956) The fitting of woven cloth to surfaces. J Text Inst 47:477-488

2. Van Der Ween F (1991) Algorithms for draping fabrics on doubly curved surfaces. Int J Numer Methods Eng 31:1414-1426

3. Long AC, Rudd CD (1994) A simulation of reinforcement
deformation during the production of preform for liquid moulding
processes. Proc Inst Mech Eng Part B J Eng Manuf 208:269-278

4. Cherouat A, Borouchaki H, Billoët JL (2005) Geometrical and mechanical draping of composite fabric. Eur J Comput Mech 14 (6-7):693-708

5. Hsiao SW, Kikuchi N (1999) Numerical analysis and optimal design of composite thermoforming process. Comp Meth Appl Mech Eng 177:1-34

6. Spencer AJM (2000) Theory of fabric-reinforced viscous fluid. Compos Part A 31:1311-1321

7. Yu WR, Pourboghrat F, Chung K, Zamploni M, Kang TJ (2002) Non-orthogonal constitutive equation for woven fabric reinforced thermoplastic composites. Compos Part A 33:1095-1105

8. Xue P, Peng X, Cao J (2003) A non-orthogonal constitutive model for characterizing woven composites. Compos Part A 34:183-193

9. Peng X, Cao J (2005) A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics. Compos Part A 36:859-874

10. Hagège B, Boisse P, Billoët J-L (2005) Finite element analyses of knitted composite reinforcement at large strain. Eur J Comput Mech 14(6-7):767-776

11. Boisse P, Gasser A, Hagege B, Billoet JL (2005) Analysis of the mechanical behaviour of woven fibrous material using virtual tests at the unit cell level. Int J Mater Sci 40:5955-5962

12. Rogers TG (1989) Rheological characterisation of anisotropic materials. Composites 20:21-27

13. Ten Thije RHW, Akkerman R, Huetink J (2007) Large deformation simulation of anisotropic material using an updated Langrangian finite element method. Comput Methods Appl Mech Eng 196(33-34):3141-3150

14. Durville D (2002) Modélisation par éléments finis des propriétés mécaniques de structures textiles: de la fibre au tissu. Eur J Comput Mech 11(2-3-4):463-477

15. Sharma SB, Sutcliffe MPF (2004) A simplified finite element model for draping of woven material. Compos Part A 35:637-643

16. Pickett AK, Creech G, de Luca P (2005) Simplified and advanced simulation methods for prediction of fabric draping. Eur J Comput Mech 14(6-7):677-691

17. Duhovic M, Bhattacharyya D (2006) Simulating the deformation mechanisms of knitted fabric composites. Compos Part A 37-11:1897-1915
18. Ben Boubaker B, Haussy B, Ganghoffer JF (2007) Discrete models of woven structures. Macroscopic approach. Compos Part B 38:498–505

19. Boisse P, Zouari B, Gasser A (2005) A mesoscopic approach for the simulation of woven fibre composite forming. Compos Sci Technol 65:429–436

20. Boisse P, Zouari B, Daniel JL (2006) Importance of in-plane shear rigidity in finite element analyses of woven fabric composite performing. Compos Part A 37:2201–2212

21. Peng XQ, Cao J, Chen J, Xue P, Lussier DS, Liu L (2004) Experimental and numerical analysis on normalization of picture frame tests for composite materials. Compos Sci Tech 64:11–21

22. Jaumann G (1911) Geschlossenes system physikalischer und chemischer differentialgesetze. Sitzungsber Akad Wiss Wien 120:385–530

23. Dafalias YF (1983) Corotational rates for kinematic hardening at large plastic deformations. J Appl Mech 50:561–565

24. Crisfield MA (1991) Non-linear finite element analysis of solids and structures, II: Advanced topics. Wiley, New York

25. Green AE, Naghdi PM (1965) A general theory of an elastic–plastic continuum. Arch Ration Mech Anal 18:251–281

26. Dienes JK (1979) On the analysis of rotation and stress rate in deforming bodies. Acta Mech 32:217–232

27. Hagege B (2004) Simulation du comportement mécanique des milieux fibreux en grandes transformations: application aux renforts tricotés. Ph.D. Thesis, ENSAM Paris

28. Yu X, Ye L, Mai Y (2005) Spurious wrinkles in forming simulations of woven fabric. Int J Form Process 8/SI:141–155

29. ten Thije RHW, Akkerman R (2008) Solutions to intra-ply shear locking in finite element analyses of fibre reinforced materials, Compos Part A (in press)

30. Boisse P, Borr M, Buet K, Cherouat A (1997) Finite element simulation of textile composite forming including the biaxial fabric behaviour. Compos Part B 28:453–464

31. Buet-Gautier K, Boisse P (2001) Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements. Exp Mech 41(3):260–269

32. Lomov SV, Verpoest I (2006) Model of shear of woven fabric and parametric description of shear resistance of glass woven reinforcements. Compos Sci Technol 66:919–933

33. Launay J, Hivet G, Duong AV, Boisse P (2008) Experimental analysis of the influence of tensions on in plane shear. Compos Sci Technol 68:506–515

34. Badel P, Vidal-Sallé E, Boisse P (2007) Computational determination of in-plane shear mechanical behaviour of textile composite reinforcements. Comput Mater Sci 40(4):439–448

35. Hamila N, Boisse P (2007) A meso-macro three node finite element for draping of textile composite preforms. Appl Compos Mater 14:235–250

36. Numisheet’93 Numerical simulation of 3-D sheet metal forming processes—verification of simulation with experiments. Edited by Makinouchi A, Nakamachi E, Onate E, Wagoner RH, Japan (1993)

37. Dumont F (2003) Experimentations et modeles de comportement de renforts de composites tisses, Ph.D. Thesis, Paris VI University