Simulations of textile composite reinforcement draping using a new semi-discrete three node finite element

Available online at www.sciencedirect.com

Composites: Part B 39 (2008) 999–1010 www.elsevier.com/locate/compositesb

Simulations of textile composite reinforcement draping using a new semi-discrete three node finite element Nahie`ne Hamila, Philippe Boisse * Laboratoire de Me´canique des Contacts et des Solides, INSA-Lyon, F-69621, France Received 10 May 2007; accepted 26 November 2007 Available online 12 January 2008

Abstract Continuous and discrete approaches for forming simulations of textile performs both involve difficulties and drawbacks. The semidiscrete method is an alternative based on a meso-macro approach. It is based on specific finite elements made of a discrete number of the components of the textile at lower scale. Their strain energy in the interpolated displacement field leads to the nodal interior loads of the element. In this paper a new three node element is proposed. It is made up of warp and weft fibres, the tensile and in plane shear energy of which are considered. The directions of the warp and weft yarns are arbitrary with regard to the element sides. That is very important in case of simultaneous multiply draping simulations and when using remeshing. The determination of the material data necessary to the simulation from standard tensile and bias tests is straightforward. A set of elementary tests and mono and multi-ply draping shows the efficiency of the approach. Ó 2008 Elsevier Ltd. All rights reserved. Keywords: A. Fabrics/textiles; C. Finite element analysis (FEA); E. Resin transfer moulding (RTM); E. Forming

1. Introduction The increase in the usage of textile composites, especially in aeronautic applications, due to their high mechanical properties leads to develop adequate design and analysis methods and especially validated simulation tools in this domain (for instance the European project ITOOL: Integrated tool for simulation of textile composites [1]). Among the different codes developed in this framework, draping simulation plays a large part. First, the simulation can give the conditions (for instance loads on the tools, initial orientation, type of material, etc.) that will make the forming of the textile reinforcement possible, and also describe possible defects after forming (wrinkles, porosities, yarn fractures, etc.). Secondly, and unique to composite forming analysis, is the need to know at any point the directions and density of the fibres after forming. These directions and densities are very important for analyzing *

Corresponding author. Tel.: +33 4 72 43 63 96. E-mail address: [email protected] (P. Boisse).

1359-8368/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesb.2007.11.008

how the composite part will behave in use, with regards to stiffness, damage, fatigue, etc. [2–4]. They also strongly influence the permeability of the textile reinforcement and must be taken into account in case of LCM processes [5,6]. In order to determine the deformed shape of draped fabrics, several codes have been developed based on geometrical approaches the so called fishnet algorithms [7–9]. 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 handmade 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 RTM process). The loads on the tools, especially on the blank-holder strongly influence the quality of the shaping operation, and therefore, need to be considered in simulations. For a physical (or mechanical) analysis of a composite forming process, the complete model must include all the

1000

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

equations for the mechanics, especially equilibrium, constitutive equations, and boundary conditions. These equations must be solved numerically, with some approximations. Finite element analysis of the composite forming process includes modelling the tools, the contact and friction between the different parts, and above all, the mechanical behaviour of the composite during forming [10]. The mechanical behaviour of textile is a multi-scale problem. The macroscopic behaviour is very much dependent on the interactions of fibres 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. 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 [11– 16]. If these models can easily be integrated in F.E. classical shell or membrane elements, the definition of the continuous constitutive model of the textile (at large strains) involves many difficulties. First there are large slidings between fibres and the continuity is questionable. Secondly the stress notion used in continuum mechanics is also questionable too because a stress vector is an elementary load on a elementary surface which, in case of fibrous material, is not well defined. Finally the update of fibre directions that play a main goal in the mechanical behaviour lead to anisotropic models that are no more orthotropic. On the other hand, the multi-scale nature of fibre reinforcements allows discrete approaches where each fibre or yarn is modelled [17–20]. Because of the large number of fibres in a fibrous reinforcement, this kind of approach is generally limited to a small domain study (woven cell, knitted loop, etc.). If the simulation of a complex part is made with such approaches, the mesoscopic modelling has to be too coarse (usually rods and springs, etc.) for an accurate result. In the semi-discrete approach [21,22], as in the discrete approach, the components at meso or micro scale are considered (yarns, fibres, woven cells etc.) and their corresponding strain energies are calculated. But in the case of the semi-discrete approach, they are part of finite elements that fix the kinematics in these discrete components as a function of nodal displacements of the elements. That leads to some continuity within the elements and especially two initially superimposed warp and weft fibers stay superimposed during the draping. That has been experimentally verified [23]. This semi-discrete approach is close to the physics by considering mesoscopic entities and accounting only for their significant mechanical properties. For instance, in the present paper, only fiber tensions and in plane shear rigidities will be taken into account. The identification of the necessary mechanical characteristics from classical tests (tensile and bias tests) is simple and straightforward. The approach proposed in this paper is based on a new three node finite element made of two fibre networks (warp

and weft), the direction of which are arbitrary with regard to the sides of the element. In the semi-discrete elements proposed in previous works, the fibres where parallel to the sides of the element [21,22]. This arbitrary orientation is essential for simultaneous multi-layer composite forming as it will be shown in Section 5.2. More generally the use of automatic meshing, remeshing and mesh refinements need this capability. The efficiency of the proposed approach will be shown in Sections 4 and 5 on elementary tests (bias test and picture frame), draping simulations and simultaneous multilayer composite forming simulations. 2. Specific mechanical behaviour of textile composite reinforcements 2.1. Tensions and in plane shear The mechanical behaviour of textiles used as reinforcement for composites is very specific [24–27]. During a forming process, there is no matrix (first stage of the RTM process for instance) or the resin is very weak and play a second order role. That is the case in draping of prepreg, where the matrix is not polymerised or in thermoforming of continuous fibre composites with thermoplastic matrix, where the matrix is liquid because heated over its fusion point. The textile composite reinforcements are made of fibres (usually carbon, glass or aramid), the diameter of which is very small (5–7 lm for carbon, 10–20 lm for aramid, 5–25 lm for glass). This diameter is extremely small in comparison to their length in case of continuous fibres. Consequently these fibres can only be submitted to tensile loads in the fibre direction. In any other direction their stiffness is close to zero. In the case of a yarn usually made of thousands of fibres, these fibres can slide relatively and the tensile stiffness is still very much larger than the others. The assembly of yarns (or fibres) in fabrics, NCF (non-crimp fabrics), knitted or braided materials leads to different mechanical properties of the resulting textile material, but the tensile stiffness in the fibre directions remain the most important. Nevertheless, for some types of analysis, it can be necessary to take the second order rigidities into account because those play a significant role. For instance, it has been shown that the in-plane shear stiffness play a main role in wrinkle development in draping simulations when the shear angle becomes large [22]. In the present work, tensile and in plane shear rigidities and corresponding strain energies will be considered. The others rigidities (bending, transverse compression etc.) will be neglected. (They could be added if necessary in some cases without modifying the main ideas developed there). 2.2. Simplified dynamic equation The textile composite reinforcement under consideration has two fibre directions (warp and weft). It can be a woven textile or a biaxial NCF (non-crimp fabric) (Fig. 1). (The

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

1001

Fig. 1. Material coordinates along the fibre directions and associated material vectors, in the cases of a woven fabric (a) and a biaxial NCF reinforcement (b).

approach could concern one direction or three and more directions of fibres but these textiles are not adapted to draping.) Material coordinates r1 and r2 along warp and weft fibre directions define the material covariant vectors k1 and k2 (Fig. 1): ox k1 ¼ 1 or

ox k2 ¼ 2 or

ð1Þ

For a fibre in direction k1, the tension on a fibre is defined by: Z T ¼ T 11 k 1  k 1 where T 11 ¼ r11 dS r ¼ r11 k 1  k 1 S1

ð2Þ r is the tensile Cauchy stress tensor, and S1 is the section of the fibre. In the same way for a fibre in direction k2: Z T ¼ T 22 k 2  k 2 where T 22 ¼ r22 dS r ¼ r22 k 2  k 2 S2

ð3Þ Considering a textile unit cell submitted to in plane shear strain such as shown Fig. 1, the warp fibres exert loads on weft fibres that can be due to weaving or to stitching in case of NCF. Their resultant is a shear torque Cs. Taking these definitions and assumptions into account, the work of the interior load in a virtual displacement field g takes the following form: nwa Z nwe Z X X p 11 p p 22 p W int ðgÞ ¼ T e11 ðgÞdl þ T e22 ðgÞdl p¼1

þ

pl

ncell X q¼1

q

p¼1

C s q cðgÞ

pl

ð4Þ

where nwa and nwe are the number of warp and weft fibres in the domain under consideration X, ncell is the number of unit cells. eaa are the axial components of the symmetrical gradient of the virtual displacement (a and b are indexes taking values 1 and 2): rs ðgÞ ¼ eab ðgÞk a  k b

ð5Þ

k1, k2 are the covariant vectors associated to k1, k2 and such as k a  k b ¼ dba . c(g) is the variation of angle between warp and weft yarns in the virtual displacement field. The form (4) of the virtual work of interior loads leads to a simplified dynamic equation that will be used in a dynamic explicit approach: nwa Z nwe Z X X p 11 p p 22 p T e11 ðgÞdl þ T e22 ðgÞdl p¼1

þ

pl

p¼1

ncell X

q

pl

C s q cðgÞ  W ext ðgÞ

q¼1

¼

Z

q€u:gdV

8g=g ¼ 0

on

Cu

ð6Þ

X

€u is the acceleration of the point P, q is the mass per volume. Cu is the part of the boundary with prescribed displacements. In this dynamic equation, the quantities are not those of continuum approaches. We have seen that those are delicate for fibrous material at large strains especially because quantities such as stresses are not well defined and because the continuity of the displacement field with the fibrous material is not clear. In Eqs. (4) and (6) the quantities such as tensions in a fibre and the torque due to the shear of a set of warp fibres on weft ones within a unit cell, are clearly defined and simple. We shall speak for this approach and for the finite element developed further of semi-discrete

1002

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

approach. The domain or element is not a continuum. It is composed of fibres in tension and shear. Nevertheless it is continuous in the sense that it is assumed there is no translation sliding between warp and weft fibres. That means that two warp and weft fibres initially superimposed at a point P (Fig. 1 for instance) will remain superimposed during the forming. That has been experimentally shown [23]. 3. Semi-discrete three node finite element 3.1. Material frames The three node triangle shown Fig. 2 is considered. In addition to the material coordinates r1, r2 along the fibres and the related material vectors k1, k2, the ‘‘natural” material coordinates (or coordinates in the reference element) n1, n2, are considered. At the three nodes of the triangle:     n1 ¼ 0  n1 ¼ 1  n1 ¼ 0    M 2 2 M 3 2 ð7Þ M 1 2 n ¼ 0 n ¼ 0 n ¼ 1 The related covariant vectors are defined: g1 ¼

ox on1

g2 ¼

ox on2

ð8Þ

The corresponding contravariant vectors are such as ga  gb ¼ dab . The material coordinates n1, n2 and r1, r2 are related. The situation of Fig. 2 is considered: k 1 ¼ AM2

k 2 ¼ BM3

M 1A ¼ a

M 1B ¼ b

ð9Þ

For a point P within the element 1

2

1

2

M 1 P ¼ n g1 þ n g 2 ¼ r k 1 þ r k 2

ð10Þ

Taking into account k 1 ¼ g1  ag2 k 2 ¼ g2  bg1

a ¼ k 1  g2 b ¼ k 2  g

1

ð11Þ ð12Þ

and the relations between the coordinates n1, n2 and r1, r2 are:        n1 1 b  r1  r1 1 b  n1 1  ð13Þ ¼  2 ¼  n 1  ab a 1  n2 a 1  r2  r2

3.2. Explicit dynamic scheme Most codes (and especially commercial ones) for material forming process simulations are based on explicit dynamic approaches [28–30] although most forming process are quasi-static. This approach will be used in the present work. It will be checked that the dynamic effects are small enough to not modify the results of the simulation. An implicit quasi static approach for fabric forming simulations has been developed in [21,23] but the explicit approach leads to a better efficiency. Within a finite element approximation, the dynamic equation in the set of degrees of freedom is: M€u þ Cu_ ¼ Fext  Fint

ð14Þ

M and C are mass and damping matrices u is the single column matrix of the degrees of freedom (three displacement components per node in the global frame). Fext and Fint are single column matrix of the components of the exterior and interior nodal loads. Eq. (14) are integrated in time using an explicit scheme using a diagonal mass matrix. Consequently, there is no need to solve a system of equations and the method is very effective. On a time step Dti, from ti to ti+1, the central difference scheme gets the solution ui+1 from ui by: uiþ1 ¼ ui þ u_ iþ1=2 Dti 1 u_ iþ1=2 ¼ u_ i1=2 þ ðDti1 þ Dti Þ€ui 2  i  €ui ¼ M1 F  Fiint  Cu_ i1=2 D ext

ð15Þ ð16Þ ð17Þ

MD is a diagonal matrix calculated from the mass matrix [31]. The time step size has to ensure the stability condition of the integration scheme, accounting for the element size and the material data [28]. The forming processes under investigations are slow and the dynamic effects are usually negligible. The explicit dynamic scheme is used because it has proved to be efficient. It is necessary to use a speed which is small enough to reduce the dynamic effects but large enough in order to reduce the computation time. This choice can be delicate. In the case of textile reinforcements, taking into account the finite element interpolation the simplified dynamic Eq. (6) can be written: gT Fext  gT Fint ¼ gT M€u þ gT Cu_ 8g=g ¼ 0 on Cu

ð18Þ

with T

g Fint ¼

nwa Z X p¼1

dl þ

p

T

11 p

e11 ðgÞdl þ

pl ncell X

nwa Z X p¼1

q

C s q cðgÞ

p

T 22 p e22 ðgÞ

pl

ð19Þ

q¼1

Fig. 2. Three nodes finite element made of textile reinforcement.

g is the single column matrix of the nodal components of the virtual displacement g. Within this explicit approach, the three node finite element presented above will be

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

defined only be by the determination of the nodal interior load Fint, or of the elementary nodal interior load Feint because Fint is assembled as:

with

Fint ¼ A Feint

W etwa int ðgÞ ¼

ð20Þ

elt

B1i1 ¼ ða  1Þk 1i

B1i2 ¼ k 1i B1i3 ¼ ak 1i ! Z n wae X p 11 T B1ij dl gij ¼ ðF etwa int Þij p¼1

A is the assembling operator of the nodal quantities from elt those on the elements to the global one. Xe is the quantity X restricted to element e. The two next sections will give the expression of the elementary nodal interior load that is, for es simplicity reason separated in tensile Fet int and shear Fint parts: ð21Þ

ð28Þ ð29Þ

pl

gij ¼ gTn Fetwa int The elementary nodal interior load are determined n wae Z X  etwa  p 11 T B1ij dl F int ij ¼

ð30Þ

pl

p¼1 es Feint ¼ Fet int þ F int

1003

3.3. Tensile elementary nodal interior load

Because of the linear interpolation, the strain interpolation terms B1ij are very simple and constants in the element. The strains and the tensions are constant in the element, consequently:

Taking the form of the virtual displacement gradient into account:



s

a

ð22Þ

etwe ¼ W etwa int ðgÞ þ W int ðgÞ

ð23Þ

gi ¼ ð1  n1  n2 Þgi 1 þ n1 gi2 þ n3 gi3

ð25Þ

gi1, gi2, gi3 are the values of gi at nodes M1, M2, M3. From Eq. (13): 1 a ½ gi1

gi2

gi3 

T

ð26Þ

Consequently Z pl

! p

T 11 B1ij dl gij



F etwe int

 ij

ð31Þ

1 ¼ kk 2 knweeB2ij T 22 2

B2i1 ¼ ðb  1Þk 2i

where i is an index taking value 1–3 and X1, X2, X3 are the components of vector X in the global frame. The interpolation functions of the three node element are the standard linear ones and the element is isoparametric. Consequently:

ogi ogi ona ¼ ¼ ½a  1 or1 ona or1

1 ¼ kk 1 knwaeB1ij T 11 2

ð32Þ

with

where nwae et nwee are the number of fibres in warp and weft directions in the element. Because the warp part W intetwa ðgÞ and weft part W intetwe ðgÞ of W et int ðgÞ are formally identical, only the calculation of the elementary nodal inteetwa rior load in warp direction Fetwa int from W int ðgÞ will be considered below.   n wae Z X p 11 p ogi ðgÞ ¼ T k W etwa dl ð24Þ 1i int or1 pl p¼1

p¼1

ij

In the same way in the weft direction:

the elementary tensile virtual work can be written:   n wae Z X og et p 11 p W int ðgÞ ¼ T  k 1 dl or1 pl p¼1   nwee Z X og p 22 p T  k 2 dl þ or2 pl p¼1

W etwa int ðgÞ ¼



b

r ðgÞ ¼ eab ðgÞk  k   og 1 og ¼  k þ  k ka  kb b a 2 ora orb

nwae X

F etwa int

ð27Þ

B2i2 ¼ bk 2i

B2i3 ¼ k 2i

ð33Þ

The elementary nodal interior load vectors due to the yarn tensions are now known. In accordance with Eq. (23), it is obtained as the sum of warp and weft contributions defined in Eqs. (31) and (32). 3.4. In plane shear nodal interior loads In the section, the in plane shear interior load vector due to in plane shear Fes int Eq. (21) is calculated in function of the shear torque resulting of actions of warp yarns on weft yarns in a unit cell and of shear strain interpolation. W es int ðgÞ ¼

ncell X

q

C s q cðgÞ ¼ gTn Fes int

ð34Þ

q¼1

The determination of the in plane shear interior load Fes int needs to give the virtual angle c(g) as a function of the nodal components of the virtual displacement gij. (The details of the calculations are given in Appendix A). The virtual shear angle c(g) (variation of angle between warp and weft directions) in the virtual displacement g is given by the gradient of g:     k1 k2 k2 k1 ð35Þ cðgÞ ¼ rg   2 þ rg   1 kk 1 k kk k kk 2 k kk k As it has been done in Eq. (26), accounting for the interpolation in the triangle, the components of the gradient of g can be expressed in function of the components of the virtual displacement gij. Consequently:

1004

cðgÞ ¼

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

 k1 k2 kk 2 k kk 1 k  2 þ B3ij þ B4ij B1ij kk 1 k kk 2 k kk 1 k kk k  2 1 k k þB2ij  g ¼ Bcij gij kk 2 k kk 1 k ij

ð36Þ

The strain interpolation terms B1ij, B2ij have been given in Eqs. (26) and (33). It is shown in Appendix A that: B3i1 ¼ ðb  1Þk 1i

B3i2 ¼ bk 1i

B3i3 ¼ k 1i

ð37Þ

B4i1 ¼ ða  1Þk 2i

B4i2 ¼ k 2i

B4i3 ¼ ak 2i

ð38Þ

The nodal interpolation of c(g) is known. Bci1 is constant in the triangle finite element, therefore:  es  F int ij ¼ ncellBcij C s ð39Þ 3.5. Material properties From the nodal interior loads defined in Eqs. (31), (32) and (39) the draping simulation can be performed using the explicit scheme Eqs. (14)–(17). As already told, one advantage of the proposed approach is that the material properties that are necessary to compute these interior loads are clearly defined. The tensions of yarns in function of their extension are obtained by a tensile test. However, in case of woven textile, there are non-linearities due to crimp and weaving and the tension in a yarn depend both on the two axial strains in warp and weft direction [32,33]. This aspect is not taken into account in the present work that is more directly applicable to NCF textiles reinforcement (Fig. 1b). Nevertheless, the crimp effect is often less important in case of woven composite reinforcement where the yarns remain few undulated. A version of this element including crimp and biaxial effect for woven reinforcement is under development. Concerning in plane shear properties, they are obtained from a picture frame test or bias extension tests; It appears

that the bias test is probably the more interesting because it avoids the spurious tensions that interfere with the picture frame results [34] and it is rather simple. In this experiment, a tensile test is performed on a rectangular specimen such as the warp and weft directions are orientated initially at ±45° to the direction of the applied tensile load (Fig. 3) [35–39]. When the specimen is stretched from L to L + d (d is the displacement of the mobile grip of the tensile machine), the possible motions between the fibres and their large tensile stiffness lead to the deformed shape shown Fig. 3c. The zone C is the active zone of the test. The warp and weft fibres have free ends. The non-sliding at crossovers and the global stretching of the specimen leads to a pure shear deformation related to d, the displacement of the mobile grip. In zone B, one yarn direction is clamped at one end, the other direction is free at both ends. The global stretching of the specimen leads to a shear strain of a value half of that of the zone C. The measure of the load on the tensile machine and of the shear angle c gives the torque Cs [40]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi  2 Dþd S B c 2 p ffiffi ffi C s ðcÞ ¼  F cD 1  Cs ð40Þ 2 2S C 2S C 2D The obtained curved is shown Fig. 4 in the case of a glass plain weave which is considered in the numerical simulations below. The characteristics of this glass plain weave (T6) are the following: yarn density = 0.25, crimp: S = 0.5%, surface density: W = 600 g/m2, linear density: w = 1200 tex (1 tex = 1 g/km). On Fig. 4, the torque is small at the beginning of the test and then it increases strongly for large shear angles. At the beginning of the relative rotation of the neighbouring yarns are only limited by the friction between warp and weft yarns. In the second part of the curve, i.e. for large shear angles, the geometry of the weaving does not allow any fur-

Fig. 3. Bias extension test (a); Initial rectangular specimen with yarns oriented at ±45° (b); Deformed specimen (c).

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

1005

Fig. 5. Simulation of the bias test. Initial and deformed meshes.

Fig. 4. Shear torque versus shear angle measured by a bias test.

Table 1 Mechanical properties of the glass plain weave Tensile stiffness in warp direction the tension is assumed to be T11 = C1e11 Tensile stiffness in weft direction The tension is assumed to be T22 = C2e22 In plane shear stiffness: C(c) = k1  c + k3  c3 + k5  c5ki are identified from Fig. 4

C1 is assumed to be constant: C1 = 46,000 N/yarn C2 is assumed to be constant: C2 = 48,000 N/yarn k1 = 0,37044 mm N k3 = 0,84105 mm N k5 = 1,03113 mm N

ther rotations and consequently the yarns are laterally compressed [33]. The angle separating the low stiff and high stiff part of the curve is called locking angle. Nevertheless it is not very precise for many fabrics and the stiffness increase can be rather progressive. This increase of shear stiffness is an important feature because it is a main cause of wrinkling [22]. The material properties of this glass fibre plain weave (shear curve shown Fig. 4) will be used in the simulations presented below. The numerical parameters are given in Table 1.

Fig. 6. Bias test: Computed and theoretical load on the tensile machine.

4. Numerical tests 4.1. Bias test The bias test presented in Section 3.5 is meshed and simulated using the three node element presented above (Fig. 5). The in plane shear behaviour Cs(c) measured Fig. 4 is used in the simulation. The computed loads on the grip of the tensile machine (F Fig. 5) is compared to the experimental load corresponding to Cs(c) i.e. given by Eq. (40). The two curves are superimposed (Fig. 6) 4.2. Picture frame

Fig. 7. Simulation of the picture frame. Initial (a) and deformed mesh (b).

The same analysis is performed in the case of the picture frame test. In this test (Fig. 7), a hinged frame with four bars with equal length is assembled in a tensile testing machine. A tensile force is applied across diagonally oppos-

ing corners of the picture frame rig causing the picture frame to move from an initially square configuration into a rhomboid. Consequently, the specimen within the picture

1006

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

Fig. 8. Picture frame: Computed and theoretical load on the tensile machine.

frame is theoretically subjected to a pure shear strain field. Fig. 7 shows the initial and deformed mesh. It can be seen on Fig. 7b that the nodal loads on the frame are well perpendicular to fibres (there is no tension in the fibres in case of pure shear). Many continuous model proposed for textile behaviour at large strain fail to this requirement. The computed loads on the grip of the tensile machine F (Fig. 7) is compared to the theoretical loads corresponding to Cs(c) [22]: ncell 2

c c 1 pffiffiffi cos  sin F ðcÞ ¼ CðcÞ ð41Þ a 2 2 2

Fig. 9. Fabric sheet draped over a circular table with stiffness: (a) woven tensile and shear stiffness, (b) null shear stiffness and (c) isotropic shear stiffness.

tions between warp and weft fibres and in plane shear is very weak. In the case of an isotropic membrane that is not possible.

The two curves are superimposed (Fig. 8). 5. Hemispherical forming 4.3. Draping on a circular table 5.1. Single ply forming The draping of a initially plane square textile sheet over a circular table is considered (The mesh is composed of 2  100  100 triangles) [18]. The proposed approach based on the three node element give a satisfactory shape after draping that involves extensive wrinkles (Fig. 9a). It can be seen that these wrinkling are well described by the membrane model. If some bending stiffness would be added, the wrinkling would be globally similar but they would be less numerous and larger as the bending stiffness increases. The draping is due to the weak shear stiffness and possible large warp weft angles variations. If the in plane shear stiffness is neglected, the draped shape is shown Fig. 9b. The draping is obtained but there is no wrinkle because the shear angle can be infinitely large as it is the case in the corners of the sheet. This solution is close to those that would be obtained by a kinematic model that appears to be unsuitable in this case. Finally Fig. 9c shows the deformed shape obtained considering an isotropic behaviour of the sheet. The draping is not obtained. It could be for instance a paper sheet that is very difficult to drape on a circular table. That show the very important role of the in plane shear stiffness in draping/forming of membrane. A textile can be shaped on a double curved surface because there are possible large rota-

Fig. 10a shows the geometry of the tools used for an hemispherical forming. The load on the blank holder is 300 N. The mesh in initial configuration is shown on Fig. 10b. It is the same for the two computed cases. In Fig. 11a and b the fibres are oriented at 0–90° and Fig. 11c and d the fibres are oriented at ±45°. The only change in the data of the simulation concerns the orientation of vectors ka. The results obtained by the two simulations and especially the shear angles (Fig. 11b and d) are in very good agreement with experimental drawings (Fig. 11a and c). 5.2. Simultaneous forming of five plies The same hemispherical forming is performed but simultaneously for five plies. The top, bottom and middle plies are oriented at 0–90° while the two intermediate plies are oriented at ±45°. The initial geometry of each ply is those of Section 5.1. The mesh is the same for each ply (Fig. 10b). That is possible with the proposed three node finite element where the directions of the fibres are arbitrary relatively to the element side. Fig. 11 shows the shear angle c in the different plies. The solution depends on the friction between

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

1007

Fig. 10. Geometry of the tools in a hemispherical deep drawing of a textile (a) and mesh used (b).

the plies which is an important point for multiply forming [41,42] (see Fig. 12). 6. Conclusions and future perspectives A semi-discrete approach for textile composite draping has been proposed. It is based on a three node finite element made up of fibres in tensions and in plane shear. The directions of the warp and weft yarns are arbitrary with regard to the element sides. Standard tensile and bias tests give directly the necessary material data. The efficiency of the element is good because only the significant quantities in textile draping are computed.

This semi-discrete approach avoids the difficulty of the determination of an equivalent continuous mechanical behaviour which in case of fibrous material at large strain has proven to be very difficult. It also very much reduce the size of the computations in comparison to discrete approaches. Nevertheless, as it has been seen in Section 3.5, the undulations due to weaving and their change in case of tension and the interaction between warp and weft yarns of woven textiles [32,33] are not taken into account in the present formulation. In practice the influence of these undulations is often small and the proposed approach can be used in these cases. Nevertherless it is better adapted

1008

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

Fig. 11. Hemispherical deep drawing of a textile reinforcement. (a) experimental, (b) computed deformed shape for 0–90° fibre orientation, (c) experimental and (d) computed deformed shape for ±45° fibre orientation.

Fig. 12. Simultaneous forming of five plies.

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

to biaxial NCF (i.e. with two layers). Taking into account the non-linearities due to weaving within a semi-discrete element is possible and a three node including this aspect is now under consideration. Bending stiffness is also neglected in the present approach. It is actually very small for most textile fabrics. Nevertheless it has been seen in Section 4 that it can modify the wrinkling shape. It will be interesting to associate the semi-discrete membrane approach developed in the present paper to bending stiffness such as taken into account in three node plate or shell elements [43–45]. The bending rigidity will have to be very specific to textile material. Then it will be possible to analyse its influence and in which case it is necessary. Acknowledgements The authors acknowledge the support of EADS aeronautics company and of the ITOOL European project. Appendix A. Interpolation of the virtual shear angle The following unit norm vector are defined: fa ¼

ka ; kk a k

fa ¼

ka kk a k

ðA1Þ

The virtual shear angle can be split in two parts: cðgÞ ¼ c1 ðgÞ  c2 ðgÞ where ca(g) is Consequently:

the

ðA2Þ angle

between

ka

and

c1 ðgÞ ¼ ðrg  f 1 Þ  f 2     og a b c1 ðgÞ ¼  kb k  k  f 1  f 2 ora   1   2 og og k k2 kk k  2 þ c1 ðgÞ ¼  k1  k1 1 2 kk 1 k kk 1 k kk k or or

ka(g). ðA3Þ ðA4Þ ðA5Þ

As seen in Section 3.4 og  k1 ¼ or1

ogi k 1i ¼ B1ij gij or1

ðA6Þ

Where B1ij are given in Eq. (28). In the same way   og og  k ¼ 2i k 1i ¼ B3ij gij 1 2 or or

ðA7Þ

With B3i1 ¼ ðb  1Þk 1i

B3i2 ¼ bk 1i

B3i3 ¼ k 1i

ðA8Þ

For the second part c2 ðgÞ ¼ ðrg  f 2 Þ  f 1     og a b c2 ðgÞ ¼   kb k  k  f 2  f 1 ora   1   1 og og k k2 kk k   c2 ðgÞ ¼   k2  k2 kk 2 k kk 1 k kk 2 k or2 or1

ðA9Þ ðA10Þ ðA11Þ

1009

And hence: og og  k 2 ¼ 2i k 2i ¼ B2ij gij or or2 Where B2ij are given in Eq. (37). In the same way   og og  k ¼ 1i k 2i ¼ B4ij gij 2 1 or or

ðA12Þ

ðA13Þ

With B4i1 ¼ ða  1Þk 2i

B4i2 ¼ k 2i

B4i3 ¼ ak 2i

ðA14Þ

References [1] ITOOL ‘Integrated Tool for Simulation of Textile Composites’, European Specific Targeted, Research Project, SIXTH FRAMEWORK PROGRAMME, Aeronautics and Space, http:// www.itool.eu. [2] Potluri P, Parlak I, Ramgulam R, Sagar TV. Analysis of tow deformations in textile preforms subjected to forming forces. Composites Science and Technology 2006;66:297–305. [3] Liu W, Drzal L, Mohanty A, Misra M. Influence of processing methods and fiber length on physical properties of kenaf fiber reinforced soy based biocomposites. Composites Part B 2007;38:352–9. [4] Mattsson D, Joffe R, Varna J. Methodology for characterization of internal structure parameters governing performance in NCF composites. Composites Part B 2007;38:44–57. [5] Bickerton S, Simacek P, Guglielmi SE, Advani SG. Investigation of draping and its effects on the mold filling process during manufacturing of a compound curved composite part. Composite Part A 1997;28:801–16. [6] Turner D, Hjelmstad K. Determining the 3D permeability of fibrous media using the Newton method. Composites Part B 2005;36:609–18. [7] Mark C, Taylor HM. The fitting of woven cloth to surfaces. Journal of Textile Institute 1956;47:477–88. [8] Van Der Ween F. Algorithms for draping fabrics on doubly curved surfaces. International Journal of Numerical Method in Engineering 1991;31:1414–26. [9] Cherouat A, Borouchaki H, Billoe¨t JL. Geometrical and mechanical draping of composite fabric. European Journal of Computational Mechanics 2005;14(6–7):693–708. [10] Boisse P. Finite element analysis of composite forming. In: Long AC, editor. Composite Forming Technologies. Woodhead Publishing; 2007. p. 46–79, Chapter 3. [11] Hsiao SW, Kikuchi N. Numerical analysis and optimal design of composite thermoforming process. Computer Methods in Applied Mechanics and Engineering 1999;177:1–34. [12] Spencer AJM. Theory of fabric-reinforced viscous fluid. Composites Part A 2000;31:1311–21. [13] Dong L, Lekakou C, Bader MG. Processing of composites: simulations of the draping of fabrics with updated material behaviour law. Journal of Composite Materials 2001;35(2):138–63. [14] Peng X, Cao J. A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics. Composites Part A 2005;36:859–74. [15] Boisse P, Gasser A, Hagege B, Billoet JL. Analysis of the mechanical behaviour of woven fibrous material using virtual tests at the unit cell level. International journal of Materials Sciences 2005;40:5955–62. [16] Akkerman R, Lamers EAD. Constitute modelling for composite forming. In: Composite Forming Technologies. Woodhead Publishing; 2007. p. 22–45, Chapter 2. [17] Durville D. Numerical simulation of entangled materials mechanical properties. Journal of Materials Science 2005;40:5941–8.

1010

N. Hamila, P. Boisse / Composites: Part B 39 (2008) 999–1010

[18] Sze KY, Liu XH. A new skeletal model for fabric drapes. International Journal of Mechanics and Materials in Design 2005;2:225–43. [19] Duhovic M, Bhattacharyya D. Simulating the deformation mechanisms of knitted fabric composites. Composites A 2006;37(11):1897–915. [20] Ben Boubaker a B, Haussy a B, Ganghoffer JF. Discrete models of woven structures. Macroscopic approach. Composites Part B 2007;38:498–505. [21] Boisse P, Borr M, Buet K, Cherouat A. Finite element simulation of textile composite forming including the biaxial fabric behaviour. Composites Part B 1997;28B:453–64. [22] Boisse P, Zouari B, Daniel JL. Importance of in-plane shear rigidity in finite element analyses of woven fabric composite preforming. Composites A 2006;37(12):2201–12. [23] Gelin JC, Cherouat A, Boisse P, Sabhi H. Manufacture of thin composite structures by the RTM process, numerical simulation of the shaping operation. Composite Science Technology 1996;56(7):711–8. [24] Scardino F. An introduction to textile structures and their behaviour, textile structural composites. In: Chou TW, Ko FK, editors. Composite Material Series. Elsevier; 1989. p. 1–26. [25] King MJ, Jearanaisilawong P, Socrate S. A continuum constitutive model for the mechanical behavior of woven fabrics. International Journal of Solids and Structures 2005;42:3867–96. [26] Long AC, Clifford MJ. Composite forming mechanisms and material characterisation. In: Composites Forming Technologies. Woodhead Publishing; 2007. p. 1, 21. [27] Lomov SV, Huysmans G, Luo Y, Parnas R, Prodromou A, Verpoest I, Phelan FR. Textile composites models: integrating strategies. Composites Part A 2001;32(10):1379–94. [28] Belytschko T. An overview of semidiscretisation and time integration procedures. In: Belytschko T, Hughes TJR, editors. Computation Methods for Transient Analysis. Elsevier Science; 1983. p. 1–65. [29] Hughes TJH, Belytschko T. A precise of developments in computational methods for transient analysis. Journal of Applied Mechanics 1983;50:1033–41. [30] Crisfield MA. Non-Linear Finite Element Analysis of Solids and Structures, 2. Wiley; 1997. [31] Zienkiewicz O, Taylor R. The Finite Element Method, 1. McGrawHill; 1989. [32] Buet-Gautier K, Boisse P. Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements. Experimental Mechanics 2001;41(3):260–9.

[33] Boisse P, Buet K, Gasser A, Launay J. Meso-macro mechanical behaviour of textile reinforcements of thin composites. Composites Science and Technology 2001;61(3):395–401. [34] Launay J, Hivet G, Duong AV, Boisse P. Experimental analysis of the influence of tensions on in plane shear. Composite Science and Technologie 2008;68:506–15. [35] Wang J, Page JR, Paton R. Experimental investigation of the draping properties of reinforcement fabrics. Composites Science and Technology 1998;58:229–37. [36] Potter K. Bias extension measurements on cross-plied unidirectional prepreg. Composites Part A 2002;33:63–73. [37] Lebrun Gilbert, Bureau Martin N, Denault Johanne. Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics. Composite Structures 2003;61:341–52. [38] Cao J, Cheng HS, Yu TX, Zhu B, Tao XM, Lomov SV, et al. A cooperative benchmark effort on testing of woven composites. In: Proceedings of the 7th international ESAFORM conference on material forming, Trondheim, Norway; 2004. p. 305–308. [39] Potluri P, Perez Ciurezu DA, Ramgulam RB. Measurement of mesoscale shear deformations for modelling textile composites. Composites Part A 2006;37:303–14. [40] De Luycker E, Boisse P, Morestin F. In: SNECMA Maia Meeting. Villaroche; 2006. [41] Gorczyca J, Sherwood J, Chen J. A friction model for use with a commingled fiberglass–polypropylene plain-weave fabric and the metal tool during thermostamping. European Journal of Computational Mechanics 2005;14(6–7):729–52. [42] Gorczyca J, Sherwood J, Chen J. A friction model for thermostamping commingled glass–polypropylene woven fabrics. Composites Part A 2007;38:393–406. [43] Belytschko T, Stolarski H, Carpenter N. A C0 triangular plate element with one-point quadrature. International Journal for Numerical Methods in Engineering 1984;20:787–802. [44] Batoz JL, Katili I. On a simple triangular Reissner–Mindlin plate element based on incompatible modes and discret constraints. International Journal for Numerical Methods in Engineering 1992;35:1603–32. [45] Boisse P, Daniel JL, Gelin JC. A C0 three-node shell element for nonlinear structural analysis. International Journal for Numerical Methods in Engineering 1994;37:2339–64.