Different approaches for woven composite reinforcement forming simulation

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.


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][2][3][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, a n dt h e r e f o r e ,n e e dt ob ec o n s i d e r e di ns i m u l a t i o n s .
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 mesoscale (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][6][7][8][9][10][11][12][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][15][16][17][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 inplane 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][8][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 )i s associated (g b Á g a ¼ d a b (a, β=1 or 2)).The (e 1 , e 2 ) basis permits an objective sum of the stress components in this frame: The stress sum (Eq. 1) is made in the (e 1 , e 2 ) basis, but the mechanical behaviour of the woven reinforcement is known in the frame of the yarns (g 1 , g 2 ).The following sets of components are considered: Assuming an orthotropic format of the constitutive matrix in the non-orthogonal frame, This relation can be written in the (e 1 , e 2 ) basis.Equation 2leads to: where T 1 ½ is the following matrix (given by Eq. 2): Consequently: D ½ is the constitutive matrix in (e 1 , 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 a and θ angles (Fig. 1).Those are calculated respectively from F (the deformation gradient tensor) and R (the polar rotation): 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][23][24]o rt h eG r e e n -Naghdi approach [24][25][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 fibre direction strictly [10].The following approach uses the rotation of the fibre denoted Δ [10,11].The rotation Δ is used to update the initial constitutive axes {k 0 } to the current constitutive axes {k t } Equation 8 leads to equations (Eq.9) that explicitly give the constitutive axes {k t } as functions of the initial constitutive axes {k 0 } and the deformation gradient F [24]: In this formulation, the fibre direction, i.e. the strong anisotropic direction, remains aligned with the first vector of {k t }.The constitutive behaviour is then fully defined at each time point.In fact, the component of the initial constitutive tensor 0 C can be computed from the traditional engineer's constants: The current constitutive tensor C can be deduced from 0 C by a rotational transport based on the fourth order rotation tensor Λ: Fig. 1 Schematic of undeformed and deformed fabric [9] Consequently C is given by: The constitutive tensor C is used in a hypo-elastic law written: D is the strain rate and s r is the objective derivative of the Cauchy stress associated to the fibre rotation Δ.The cumulated tensorial strain tensor ɛ and stress tensor σ associated with such an objective derivative are: 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 Δ for the Fig. 3 Simulation result of a bias extension test on a NCF using a mesoscopic or discrete approach [16] 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] 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.

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.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: where ɛ h ðÞ¼r s h ¼ " αβ h ðÞ h α h β is the symmetrical gradient in the virtual displacement η (a and β are indexes taking value 1 or 2).h 1 , h 2 are the contravariant vectors related to the unit vector in warp and weft directions h 1 and h 2 .L 1 and L 2 are the lengths of the warp and weft yarns in the mid-plane of the fabric.g(η) 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, p Q, means that the quantity Q is considered for the woven unit cell number p. ρ is the mass per volume of the fabric Ω. W ext (η) is the virtual work of the exterior prescribed loads.Γ 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 (" 11 , " 22 ), T 22 (" 11 , " 22 ) and C(g).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 g) 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(g)i s made by default.
The tension surfaces (T 11 (" 11 , " 22 ) and T 22 (" 11 , " 22 )) and of the shear curve (C(g)) 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 1 Mechanical properties of the glass plain weave ("Deep drawing of a square box.Influence of the in plane shear stiffness" section) Tensile stiffness in warp direction the tension is assumed to be T 11 =C 1 " 11 C 1 is assumed to be constant: C 1 =46,000 N/yarn Tensile stiffness in weft direction the tension is assumed to be T 22 =C 2 " 22 C 2 is assumed to be constant: C 2 =48,000 N/yarn In plane shear stiffness: Eq. 5 the nodal interior load components F e int ÀÁ s are obtained from the strain interpolation matrix: 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  2 Mechanical properties of the unbalanced fabric ("Hemispherical forming of an unbalanced composite fabric" section) Tensile stiffness in warp direction the tension is assumed to be T 11 =C 1 ɛ 11 C 1 is assumed to be constant: C 1 =50 N/yarn Tensile stiffness in weft direction the tension is assumed to be T 22 =C 2 ɛ 22 C 2 is assumed to be constant: In plane shear stiffness: .03 mm N/rad, k 2 =0.095 mm N/rad and, g c =0.5 rad (29°) 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 asymmetrical.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 weft /l warp 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 semidiscrete approaches the directions of the yarns are "naturally" 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.

Fig. 4
Fig. 4 Four node finite element made of woven cells

Fig. 8
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 25mm N/rad, k 2 =4.08 mm N/rad and, g c =0.72 rad (41°)