- Email: [email protected]
Hyperelastic modelling for mesoscopic analyses of composite reinforcements
Composites Science and Technology 71 (2011) 1623–1631
Contents lists available at ScienceDirect
Composites Science and Technology journal homepage: www.elsevier.com/locate/compscitech
Hyperelastic modelling for mesoscopic analyses of composite reinforcements A. Charmetant, E. Vidal-Sallé, P. Boisse ⇑ Université de Lyon, LaMCoS, INSA-Lyon, F-69621, France
a r t i c l e
i n f o
Article history: Received 21 April 2011 Received in revised form 24 June 2011 Accepted 3 July 2011 Available online 20 July 2011 Keywords: A. Fabrics/textiles B. Non-linear behaviour C. Anisotropy C. Finite element analysis Hyperelasticity
a b s t r a c t A hyperelastic constitutive law is proposed to describe the mechanical behaviour of fibre bundles of woven composite reinforcements. The objective of this model is to compute the 3D geometry of the deformed woven unit cell. This geometry is important for permeability calculations and for the mechanical behaviour of the composite into service. The finite element models of a woven unit cell can also be used as virtual mechanical tests. The highlight of four deformation modes of the fibre bundle leads to definition of a strain energy potential from four specific invariants. The parameters of the hyperelastic constitutive law are identified in the case of a glass plain weave reinforcement thanks to uniaxial and equibiaxial tensile tests on the fibre bundle and on the whole reinforcement. This constitutive law is then validated in comparison to biaxial tension and in-plane shear tests. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Their good mass-mechanical properties ratio and their good corrosion and fatigue properties make textile composites essential materials when saving mass is a main issue. These materials can be manufactured by Liquid Composite Molding (LCM) processes [1– 3]: the textile reinforcement is preformed in a first stage, then the resin is injected in the preform. The reinforcement is made of continuous fibres grouped into yarns (3000–48,000 fibres in a carbon yarn, 1000–12,000 fibres in a glass yarn). These yarns are combined following textile architectures: woven, knitted braided or non-crimped fabrics. These assemblies give the textile reinforcements a multi-scale nature (Fig. 1). This paper focuses on woven fabrics, but the model which will be proposed for the yarn would be similar for other architectures. The macroscopic scale refers to the whole component level, with dimensions ranging from ten centimetres to several metres. At this level, a woven fabric can be seen as a continuous material whose behaviour is very specific and includes high anisotropy. At the mesoscopic scale, the woven reinforcement is seen as interlaced tows, respectively the warp and the weft (or fill) yarns. Consequently, the working scale corresponds to the yarn dimension, typically one to several millimetres. For periodic materials, mesoscopic models consider the smallest elementary pattern, which can represent the whole fabric by translations. That domain is called the representative unit cell (RUC) [4–6]. Finally, the rein-
⇑ Corresponding author. Tel.: +33 4 72 43 63 96; fax: +33 4 72 43 85 25. E-mail address: [email protected] (P. Boisse). 0266-3538/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compscitech.2011.07.004
forcement can be analysed at the microscopic scale. In this case each fibre is considered as a 3D beam interacting with its numerous neighbours [7,8]. At the microscopic level, the characteristic dimension is about one to several micrometres. This is the only scale at which the material is actually continuous. During the forming process, the woven unit cell reaches its position and geometry in the final composite. This geometry is important regarding the resin flow [9,10] and the final mechanical properties of the composite part. The changes of geometry of the fibre bundles are analysed using mesoscopic scale modelling. The studied domain is then the representative unit cell, which gives the solid skeleton for permeability evaluation. In addition to permeability calculation, mesoscopic analyses can be used as virtual mechanical tests: they allow the mechanical properties of a woven reinforcement to be obtained without performing the experiments and hence without necessarily manufacturing the reinforcement under consideration [4]. In mesoscopic computations the yarn is considered as a continuum but with a very specific behaviour which accounts for his fibrous and discontinuous composition. Some mesoscopic analyses of the woven unit cell submitted to tension, in-plane shear and compression have been set up based on a hypoelastic behaviour of the yarn [4,5,11,12]. These hypoelastic models require the use of an objective derivative, i.e. of a derivative in a frame as close as possible of the fibre orientation. Despite being fairly frequently employed, these models have drawbacks, especially they do not allow a complete recovering after a closed loop loading path. The objective of the development of a hyperelastic constitutive model is to avoid these drawbacks. In addition, a hyperelastic model has been developed recently to analyse the mechanical behaviour of macroscopic textile composite reinforcements modelled as membranes with
1624
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Macroscopic scale Mesoscopic scale Microscopic scale Fig. 1. Different scales for a woven composite reinforcement.
satisfactory results [13]. Nevertheless the model necessary to analyse the deformation of the fibre bundle is quite different since it must be 3D with a single fibre direction while the previous macroscopic model is 2D with two fibre directions. Furthermore the transverse deformation of the yarn (i.e. its compression and shape change), which is essential to the yarn behaviour, cannot be described by a macroscopic membrane approach. An important research effort has been devoted to the development of hyperelastic models through the last decades. Isotropic models have been developed first, to model the behaviour of materials which can be subjected to very large elastic deformation, such as rubber–like materials. These models are usually sorted into two classes, depending on whether the strain energy depends on invariants of the strain tensor in use (e.g. Mooney–Rivlin models [14]), or on principal elongations (e.g. Ogden models [15]). The development of hyperelastic constitutive models for anisotropic materials appeared more recently. They are used to model the mechanical behaviour of strongly anisotropic materials which can undergo large elastic deformations, e.g. reinforced elastomers [16–18] and some biomaterials (biological tissues [19], tendons and ligaments [20]). The formulation of anisotropic hyperelastic models is usually based on structural tensors [21–23] which describe the local orientation (i.e. the preferred directions) of the material. In this paper an anisotropic hyperelastic model is proposed and validated for the mechanical behaviour of textile composite reinforcements. It is based on specific physically motivated invariants. First the specificities of the yarn behaviour are analysed in order to formulate the different assumptions of this work. The mechanical behaviour model of the yarn is then presented, and the identification of its parameters is performed in the case of a glass plain weave, by use of uniaxial and biaxial tensile tests. Finally the model is used to simulate the in-plane shear behaviour of the woven unit cell, and the obtained results are compared to experimental data. 2. Mechanical behaviour of the fibre bundle and assumptions 2.1. Transverse isotropy The yarn is an assembly of fibres oriented approximately in the same direction (Fig. 2). The fibres are assumed to be numerous and compact enough in order to avoid independent movements. In this case, it is possible to consider the yarn as a continuum, and consequently to define a continuous material equivalent to the fibre bundle. In a recent study [24], an X-ray tomography technique was combined to a mechanical test on a fibre bundle. This ap-
Fig. 2. Transverse isotropy of the yarn in a textile reinforcement. X-ray tomography imaging.
proach allowed an accurate observation of the motion of fibres inside the yarn, and showed that such a continuous approach is wellfounded. The homogenized material has a preferred direction which is the fibre one. The spatial distribution of fibres inside the yarn has been analysed on deformed cross sections and, even if this observation may differ for other types of fabrics [25], it has been concluded that this distribution is isotropic [5] for the considered fabric. The ‘‘yarn’’ material will then be assumed to be transversely isotropic.
2.2. Deformation modes A phenomenological approach based on experimental observations is used to define the constitutive law. It is consequently necessary to describe the specific deformation modes of the yarn, which are related to the behaviour of the fibres and to their interactions at the microscopic scale. Four deformation modes are considered: the elongation in the fibre direction, the compaction in the transverse section of the yarn, the distortion (shear) in the transverse section and the shear along the fibre direction (Fig. 3). The last one (longitudinal shear) is a measure of the angle between the covariant cross section and the initial cross section. It is then compound of both longitudinal shear deformations. This deformation mode mainly controls the bending rigidity of the yarn. The experimental analysis of the tensile behaviour of the yarn is quite easy, so its identification can be performed directly from a tensile test on the yarn. This is not the case for the three other deformation modes. The weaving gives cohesion to the fibre bundle in the yarn and it is not easy to define elementary tests for the compaction and shear on a single yarn because of the lack of this cohesion if the yarn is extracted from the fabric. An inverse approach is used in the present work for the determination of the material parameters in transverse compression and shear (modes
1625
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
(a)
(b)
(c)
(d)
Fig. 3. Deformation modes of the yarn: (a) elongation, (b) compaction of the cross section, (c) distortion of the cross section and (d) longitudinal shear. The arrow indicates the direction of fibres.
(b–d) Fig. 3). It is performed on uniaxial and equibiaxial tensile tests because these three modes are present in these cases. 2.3. Uncoupling assumptions The proposed model does not take into account the probably existing couplings between the behaviours corresponding to the deformation modes defined above. It is then assumed that, in the range of deformation simulated and for the considered loading cases, the strain energy function can be split into four parts, each being a function of only one deformation mode. 3. Hyperelastic constitutive equation 3.1. Finite displacement kinematics
~ CQ T Þ; ~ ¼ wðQ wðCÞ
8Q 2 G
ð6Þ
The symmetry group of an isotropic material includes all the rotations and reflections (i.e. G ¼ O3 , the group of isometries of R3 ). Con~ is mathematically isotropic. For an sequently, the strain energy w ~ is mathematically anisotropic. In anisotropic material, GˆO3 so w such a case, it has been shown [21–23] that an isotropic function ^ exists which depends on the right Cauchy–Green strain tensor w C and on N structural tensors Gi defining the symmetry group, such as:
^ ~ ¼ wðC; G1 ; :::; GN Þ wðCÞ
Let X be the initial position of a material point and x its position in the current configuration. v is the transformation from the initial to the current configuration: x ¼ vðXÞ. The deformation gradient tensor is then:
ð7Þ
The representation theorems for transversely isotropic functions of scalars, vectors and/or symmetric tensor variables define a finite ~ such as: number of scalar functions ðI1 ; . . . ; In Þ and a function w _
^ wðC; G1 ; . . . ; GN Þ ¼ wðI1 ; . . . ; In Þ
@x F¼ @X
ð1Þ
This tensor allows the transformation of any vector A in the initial configuration into a vector a of the current (deformed) configuration. The Cauchy–Green deformation tensor C ¼ F T F, from which the variation of the length of A can be deduced, is defined as follows:
kak2 ¼ ðF AÞ ðF AÞ ¼ A ðF T FÞ A ¼ A C A
ð2Þ
The invariants of tensor C are usually defined by:
I1 ¼ traceðCÞ;
1 I2 ¼ ðtraceðCÞ2 traceðC 2 ÞÞ; 2
I3 ¼ detðCÞ
ð3Þ
3.2. Hyperelastic equations
ð4Þ
where F is the deformation gradient tensor and S is the second Piola–Kirchhoff stress tensor. To ensure that this constitutive equation fulfils the principle of ~ of the right Cauchy Green material frame indifference, a function w strain tensor C must exist as [26]:
~ ~ T FÞ ¼ wðCÞ ¼ wðFÞ wðF
Number n depends on the number and type of arguments of the isotropic function under consideration. Finally, the strain energy function will be defined in the form (8). For sake of simplicity, it will be denoted w in the sequel. The symmetry group of a transversely isotropic material is characterised, in the initial configuration, by an unit vector M which represents the preferred direction. This allows the definition of a structural tensor M ¼ M M. Then, the strain energy potential becomes:
ð9Þ
where I1, I2, I3 are the standard invariants of C defined in (3). II4 and II5 are mixed invariants defined from the structural tensor M by:
II4 ¼ C : M
and II5 ¼ C 2 : M
ð10Þ
The second Piola–Kirchhoff stress tensor is then defined as:
In order to define a hyperelastic material, an elastic strain energy potential per unit volume w must exist which only depends on the current strain state [26]. As the behaviour is reversible, it does not dissipate energy and naturally fulfils the Clausius–Duhem inequality. These assumptions allow the constitutive equation of a hyperelastic material to be written as:
@wðFÞ @C
ð8Þ
w ¼ wðI1 ; I2 ; I3 ; II4 ; II5 Þ
The local change of volume during deformation is described by the pffiffiffiffi Jacobian J ¼ I3 ¼ detðFÞ.
SðFÞ ¼ 2
Consequently, the strain energy potential can be defined directly, ~ using wðCÞ. Moreover, it must satisfy the principle of material symmetry, i.e. the strain energy must be invariant under any transformation belonging to the symmetry group G of the material:
ð5Þ
S¼2
@w @C
@w @I1 @w @I2 @w @I3 @w @II4 @w @II5 þ þ þ þ ¼2 @I1 @C @I2 @C @I3 @C @II4 @C @II5 @C
! ð11Þ
4. Constitutive equation for a yarn In this section, a hyperelastic model is proposed which describes the behaviour of the yarn in textile composite reinforcements. First, ‘‘physically based’’ invariants are defined which allow the description of the different deformation modes of the yarn introduced in Section 2. Then, each mode will be associated a strain energy function from which the stresses in the material can be derived.
1626
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
4.1. Invariants and their corresponding deformation modes An approach based on a decomposition of the deformation gradient [27] is used to determine physically based invariants of the transformation for each deformation mode of the fibre bundle. An orthonormal basis B ¼ fM; N 1 ; N 2 g is defined which allows writing the components of the deformation gradient F of transversely isotropic materials into the form:
2
fm 6 ½FB ¼ 4 0 0
fm1
fm2
3
0
0
0
32
1
½F elong B
2
½F comp B
1
0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 f11 =f22 40
32
0
1
fm1 fm
fm2 fm
3
76 7 0 5 40 1 0 5 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 0 f22 =f11 0 0 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} ½F dist B
ð13Þ
½F sh B
where F elong , F comp , F dist and F sh respectively describe the elongation of the yarn along the fibres direction, the compaction and distortion of the yarn in the transverse direction, and the longitudinal shear of the yarn. Each part of the decomposition can be characterized by a single quantity, except F sh which requires two quantities. Finally, five parameters are needed to describe the deformation with this approach:
aelong ¼ fm ; acomp ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f2 þ f2 ash ¼ m1 2 m2 ; fm
pffiffiffiffiffiffiffiffiffiffiffi f11 f22 ;
tanðcÞ ¼
adist ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi f11 =f22 ;
fm1 fm2
ð14Þ
This is consistent with the fact that five invariants are needed to write the strain energy potential (9) of hyperelastic transversely isotropic materials. With these notations, the components of the transformation gradient can be written as:
2
sffiffiffiffiffiffiffiffiffiffiffiffiffi I5 1 I24
ð17Þ
These invariants are used in the sequel to define the strain energy functions associated to each deformation mode:
ð18Þ
4.2. Strain energy functions associated to each deformation mode
3 0 0 p ffiffiffiffiffiffiffiffiffiffiffi 6 76 0 7 f11 f22 ½FB ¼ 4 0 1 0 5 4 0 5 pffiffiffiffiffiffiffiffiffiffiffi 0 0 1 0 0 f11 f22 |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} fm
Ish ¼
ð12Þ
As previously mentioned, M is a unit vector which direction is the preferred direction of the material (i.e. the direction of fibres) in the initial configuration. A multiplicative split of the matrix (12) can be made in order to separate the parts corresponding to each deformation mode of the yarn:
2
Idist
wðCÞ ¼ wðIelong ; Icomp ; Idist ; Ish Þ
7 0 5 f22
f11
1 1 I3 ; lnðI4 Þ; Icomp ¼ ln 2 4 I4 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 2 1 I 1 I 4 I5 I1 I 4 I 5 pffiffiffiffiffiffiffiffi ¼ ln @ pffiffiffiffiffiffiffiffi þ 1A; 2 2 I3 I4 2 I3 I4
Ielong ¼
3
aelong ash aelong cosðcÞ ash aelong sinðcÞ 7 acomp adist 0 5 0 acomp =adist
6 ½FB ¼ 4 0 0
sffiffiffiffiffiffiffiffiffiffiffiffiffi I ash ¼ 52 1 I4
nl w elong ðIelong Þ if welong ðIelong Þ ¼ wlin ðIelong Þ if elong
Ielong 6 I0elong
ð19Þ
Ielong > I0elong
The strain energy functions corresponding to the non-linear and linear parts of the elongation behaviour of the yarn are respectively defined by the following equations: 2 3 wnl elong ðIelong Þ ¼ a0 þ a1 I elong þ a2 Ielong þ a3 Ielong
ð20Þ
2 wlin elong ðIelong Þ ¼ b0 þ b1 I elong þ b2 I elong
ð21Þ
Introducing the initial yarn section S0 (assumed to be constant along the yarn), its initial stiffnesses K 0elong , its stiffness in the linear part Kelong, and the continuity condition between the two functions lin wnl elong et welong , leads to:
wnl elong ¼
wlin elong ¼
ð15Þ
K 0elong 2 K elong K 0elong 3 Ielong þ Ielong 2S0 6S0 I0elong
ð22Þ
K elong K 0elong 0 2 K elong K 0elong 0 Ielong Ielong Ielong 6S0 2S0 K elong 2 þ I 2S0 elong
ð23Þ
Using (11), the second Piola–Kirchhoff stress tensor is given by:
The quantity c is a measure of the coupling between the two shear modes, namely distortion and longitudinal shear [27,28]. Under the uncoupling assumptions formulated before, the behaviour of the yarn is then assumed not to depend on this quantity. The standard invariants (9) can then be expressed as functions of the quantities (14). Inverting the obtained expressions leads to:
vffiffiffiffiffiffiffiffi ffi usffiffiffiffi u I pffiffiffiffi aelong ¼ I4 ; acomp ¼ t 3 ; I4 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 uI I I I1 I4 I5 t14 5 p pffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffi ffi þ adist ¼ 1; 2 I3 I4 2 I3 I4
4.2.1. Elongation strain energy A tensile test on a single yarn allows the identification of its tensile behaviour. Such a test highlights a response of the yarn in tension which can be close to a linear behaviour (see Fig. 4). Nevertheless, the slight misalignment of the fibres inside the tow can induce an initial non linearity: as the fibres do not have exactly the same length, the overall tension of the yarn corresponds to a non simultaneous tension of the fibres, and consequently to a progressively stiffening behaviour of the yarn [29]. In order to describe both linear and non-linear parts of the response of the yarn in tension, the corresponding strain energy potential is defined piecewise:
ð16Þ
Amongst the previous quantities, only ash vanishes when the material is not deformed (i.e. when F ¼ I). After normalizing the other quantities, the following set of invariants is defined:
Selong
0 K elong K 0elong 2 K elong 1 S0 Ielong þ 2S0 I0elong Ielong ¼ M I4 K elong K 0elong 0 K I þ elong Ielong 2S0
elong
S0
if
Ielong 6 I0elong
if
Ielong > I0elong
ð24Þ
Four parameters have to be identified: S0, I0elong , K 0elong and Kelong. The first one will be determined by measuring the geometrical characteristics of the fabric. To avoid stiffness errors, S0 must exactly match the initial cross-sectional area defined in the finite element model. The three remaining parameters will be identified using the result of the tensile test. 4.2.2. Compaction strain energy The strain energy associated to the transverse compaction of the yarn is linked to several physical phenomena at the microscopic scale (mostly bending and rearrangement of fibres) which
1627
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Identified parameters for elongation
200
S0 (mm²) I elong 0 0 (N) K elong K elong (N)
150
Simulation
Effort (N)
1.0496 0.00134 11660 36714
Experience
100 50 0 0
0.002
0.004
0.006
Deformation of the yarn Fig. 4. Identification of yarn elongation behaviour.
are difficult to study separately. As in [4,30–32], a power-based strain energy function is then proposed:
p K comp jIcomp j wcomp ðIcomp Þ ¼ 0
if
Icomp 6 0
if
Icomp > 0
ð25Þ
ð26Þ
The parameters to be identified are then Kcomp and p. As a direct experimental characterization of this behaviour is difficult, an inverse method is set up, based on equibiaxial tensile tests, to identify both parameters. An equibiaxial tensile test is used rather than a crushing test because the later has been found to be very sensitive to the frictional coefficient between the crushing tools and the fabric, which may make the identification less reliable. Moreover, crushing simulations computation time is much larger, which significantly increases the duration of the identification procedure. 4.2.3. Distortion strain energy As well as those of the compaction behaviour, the parameters of the behaviour of the distortion in the transverse section can not be determined by a direct experiment. In the present work, it has been assumed that the distortion stiffness is constant which leads to the following strain energy function:
wdist ðIdist Þ ¼
1 K dist I2dist 2
ð27Þ
Constant Kdist will be identified at the same time as the parameters of Eq. (25), using the same inverse method. The second Piola–Kirchhoff stress tensor associated to this distortion strain energy is then given by:
Sdist ¼ 2K dist Idist
2I4 I ðI1 I4 I5 ÞC 1 þ I1 þ II54 M 2ðC M þ M CÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 ðI1 I4 I5 Þ2 4I3 I4
1 K sh I2sh 2
ð29Þ
The second Piola–Kirchhoff stress tensor is then given by:
It is assumed that no energy is required for yarn section expansion. It is also assumed that during the forming of a textile composite reinforcement, the yarn compaction stiffness never reaches the bulk stiffness of glass (i.e. the voids in the yarn do not collapse altogether). The second Piola–Kirchhoff stress tensor, corresponding to the yarn compaction, derives from Eq. (25):
p 1 Scomp ðIcomp 6 0Þ ¼ K comp jIcomp jp1 C 1 M 2 I4
wsh ðIsh Þ ¼
ð28Þ
4.2.4. Longitudinal shear strain energy The longitudinal shear strain energy of the yarn is linked to the sliding of its constitutive fibres in their preferred direction. A linear elastic behaviour will be assumed for this deformation mode, which induces the following strain energy potential:
" # 1 1 2I5 Ssh ¼ K sh 2 ðC M þ M CÞ 3 M 2 I4 I4
ð30Þ
During a tensile test on a textile composite reinforcement, this deformation mode is activated: during the non-linear part of the tensile curve that corresponds to a change of undulation of the yarns, before reaching the maximum tensile stiffness, longitudinal shear and elongation are the only working modes, whereas in the linear part, only the longitudinal elongation works (Fig. 7). The identification of Ksh is then possible using an inverse method on the uniaxial tensile test as soon as the longitudinal elongation behaviour is known. 4.3. Combination of all contributions Finally, the whole constitutive equation is the summation of all contributions discussed before:
S¼2
@welong @Ielong @wcomp @Icomp @wdist @Idist @wsh @Ish þ þ þ @Ielong @C @Icomp @C @Idist @C @Ish @C
!
¼ Selong þ Scomp þ Sdist þ Ssh ð31Þ Table 1 summarizes the parameters which have to be identified for each deformation mode, as well as the experimental test used for their identification. The proposed constitutive equation has been implemented as a user material subroutine VUMAT in the ABAQUS/Explicit finite element code. 4.4. Studied reinforcement and finite element model of the yarn For sake of simplicity, the parameter identification has been realized on a glass plain weave which main geometrical parameters are summarized in Table 2. The plain weave presents the great advantage to have a rather simple unit cell: it is possible to obtain a realistic geometrical model of the yarn, using only straight lines and circular arches, avoiding yarn interpenetrations [33]. Fig. 5 shows the geometry used in the present work. The elementary unit cell has several symmetry axes and planes. When the boundary conditions imposed to the unit cell are compatible with those symmetries, it becomes possible to reduce the pattern to be studied up to 1/8 of the unit cell (Fig. 6). For example, it is the case of biaxial tension and crushing. Nevertheless the whole unit cell has to be used in order to simulate in-plane shear. Contact and frictional sliding between the tows is treated using Coulomb friction with a coefficient of 0.3 identified in [34].
1628
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Table 1 Summary of the parameters and their identification approaches. Deformation mode
Elongation
Compaction
Number of parameters Identification mode
4 Geometric observations, tensile test on the yarn
2 1 Equibiaxial tensile test, iterative approach
Plain weave
Yarn width (mm)
Warp: 3.2 Weft: 3.1
Yarn density (yarn/mm)
Warp: 0.251 Weft: 0.248
Crimp (%)
Warp: 0.5 Weft: 0.54
Surface weight (g/m2)
600
-3
Weave type
Energies (mJ.mm )
Table 2 Geometrical characteristics of the studied glass plain weave.
Distortion
Longitudinal shear 1 Tensile test on the fabric
1.00E+00 1.00E-01 1.00E-02 1.00E-03 1.00E-04 1.00E-05 1.00E-06 1.00E-07 1.00E-08 1.00E-09 1.00E-10 1.00E-11 1.00E-12
Welong Wsh Wcomp Wdist
0
0.02
0.04
0.06
0.08
0.1
0.12
Deformation of the fabric Fig. 7. Evolution of strain energies associated to each deformation mode during an uniaxial tensile test on the reinforcement.
(a)
(b)
Circle arches
Straight segments
Circle arches
Fig. 5. Real geometry obtained using X-ray tomography imaging: (a) and corresponding model and (b) for a unit cell of plain weave.
yarn bending). Once the parameters of the elongation behaviour of the yarn have been identified, it becomes, then, possible to identify the parameter Ksh of Eqs. (29) and (30). This method allows determining the value Ksh = 3 MPa. On Fig. 7, the strain energies associated to each deformation mode during a uniaxial tensile test on the fabric are plotted. Since the energies of elongation and longitudinal shear are much greater than those associated to compaction and distortion, this method of identification using this kind of test seems well-founded.
5. Parameters identification
5.3. Identification of compaction and distortion behaviours
5.1. Identification of the behaviour in elongation
As mentioned before, the parameters Kcomp, Kdist and p from Eqs. (25) and (27) are then identified using an inverse method: they are optimised so that the simulated behaviour of the unit cell matches its experimentally observed behaviour. To ensure an accurate identification, the experimental test of reference must mainly depend on these deformation modes of the yarn. In this paper the experimental test of reference chosen is equibiaxial tension which, as seen afterwards, satisfies the former condition. An equibiaxial tension test consists in imposing the same elongation to both warp and weft directions of the woven reinforcement. Such a test induces a significant amount of compaction and distortion: because the tensile stiffness of the yarn is very high, a fabric which undergoes equibiaxial tension starts adapting itself to the deformation by reducing its crimp (i.e. the undulations of the yarns). This crimp reduction is made possible by both compac-
The initial section of the yarn, S0, is given by mesoscale observations of the fabric, using X-ray tomography for instance. The three other parameters, namely I0elong , K 0elong and Kelong, are identified thanks to a tensile test on a single yarn. As the yarn may have been damaged during the weaving, this test must be performed on a yarn extracted from the fabric rather than a yarn before weaving. The parameters obtained by this method for the glass fibre yarn, as well as the experimental and identified curves, are shown on Fig. 4. 5.2. Identification of the longitudinal shear behaviour During an uniaxial tensile test on the fabric, the main deformation modes of the yarn are elongation and longitudinal shear (i.e.
(a)
(b)
(c)
Fig. 6. (a) Plain weave. (b) Elementary unit cell for in-plane shear. (c) 1/8 of the unit cell for biaxial tension.
1629
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Identified parameters for compaction
140
K comp (N) pcomp
120
Experience (warp direction)
Identified parameter for distortion K dist 0.06
Effort per yarn (N)
17800 7.72
Experience (weft direction) Initial simulation (10000,5,0.05)
100
Loop 1 (9300,6.82,0.0542) Loop 2 (10300,7.16,0.0561)
80
Loop 3 (17800,7.72,0.0608)
60 40 20 0 0
0.001 0.002 0.003 0.004 0.005 0.006 Deformation of the fabric
Fig. 8. Optimisation of the parameters of compaction and distortion using equibiaxial tension.
200 180
Effort per yarn (N)
tion and distortion. That is why the non-linear part of the curve of equibiaxial tension (Fig. 8) is representative of the behaviour of the yarn in its transverse plane. The inverse method which has been set up consists in minimizing a least square estimation of the error between the experimental curve of equibiaxial tension, and the curve obtained by simulation. An iterative Levenberg–Marquardt algorithm [35,36] has been used to determine the three parameters of compaction and distortion. Fig. 8 illustrates the successive curves obtained during this optimisation. The proposed method allows to obtain a result of simulation very close to the observed behaviour, with only few iterations. Nevertheless, the existence of local minima of the error estimation induces that the initial parameters must be close enough to the final solution so that the algorithm converges. This imposes a first trial and error approach before starting the optimisation procedure.
Experience Simulation
160
Yarn k=0.5 free
140 120
k=1
100 80
k=2
60 40 20 0 0
0.002
0.004
0.006
0.008
Deformation of the fabric Fig. 9. Biaxial tension curves obtained after identification of the parameters.
6. Validation of the hyperelastic law In this section, the previous hyperelastic constitutive law and the identification of its parameters are validated: the parameters previously identified are used to perform simulations in other cases and the results are compared to experimental values. A first validation is established by comparing simulation and experimental results of biaxial tension with different warp/weft elongation ratios. Fig. 9 presents the whole network of experimentally and numerically obtained curves. Remind that the curves ‘‘yarn’’, ‘‘free’’ and ‘‘k = 1’’ have been used for identification. The validation can then be performed by comparing simulation and experimental curves for the cases ‘‘k = 0.5’’ and ‘‘k = 2’’. The two numerically obtained curves show good agreement as regard as experimental curves, which allows a first validation of the proposed hyperelastic constitutive law. Biaxial tension tests with k = 0.5 and k = 2 elongation ratios are, kinematically, relatively close to the k = 1 test used for the identification. That is why another validation is performed using an experimental test which kinematics is very different from biaxial tension. Then, a simulation of the shear behaviour of the unit cell has been performed and compared to the results of the picture frame test [37–40]. In order to compare the meso-scale simulation to the macro-scale picture frame test, it is assumed that the fabric and the unit cell have the same mechanical behaviour, i.e. that the shear deformation in the fabric during a picture frame test is uniformly distributed. The periodicity and boundary conditions imposed to the unit cell in this simulation of shear are similar to
those described in [5]. The initial and deformed configurations of the unit cell are shown on Fig. 10. Two criteria are used to compare the results: one which allows validating the overall efforts induced in the unit cell and one to validate its deformed geometry. In order to compare the results obtained by different measurement devices, or to compare the inplane shear behaviour of different types of fabrics, a quantity must be defined which characterize the efforts in the unit cell. Different quantities have been introduced in previous papers. In the sequel the torque Cs per unit initial area introduced in [40] is used:
Cs ¼ F
cosðcÞ L2fabric 2cosðp=4 c=2Þ
Lframe
ð32Þ
with F the effort measured on the frame, Lframe the side of the frame, Lfabric the side of the square of fabric inside the frame and c the angle variation between warp and weft directions (i.e. the shear angle). On the other hand, an energy-based approach [40] allows the calculation of this torque per unit initial area (32) from the simulation:
Cs ¼
_ W Su c_
ð33Þ
_ is the time-derivative of the external work, Su the initial where W surface of the unit cell and c_ the time-derivative of the shear angle. The comparison of the torque obtained from experiments and simulations is presented on Fig. 11a. On this figure, the shear torques computed during the simulation and measured during the experi-
1630
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
0.7 Picture frame test
0.6
Simulation
0.5 0.4 0.3 0.2 0.1 0 0
(a)
20
40
60
Shear angle (deg)
Relative broadness of the yarns
Torque per unit initial area (N)
Fig. 10. Unit cell shear simulation: (a) undeformed configuration and (b) deformed configuration.
1 0.8 0.6 0.4
Experiment
0.2
Simulation
0 0
10
(b)
20
30
40
50
60
Shear angle (deg)
Fig. 11. (a) Evolution of normalized shear force during numerical and experimental tests. (b) Evolution of yarn broadness during numerical and experimental tests.
ment are close. During the preparation of shear samples for experiments, a slight tension of the fabric is required to ensure a good positioning of the fabric in the picture frame. This tension induces stronger interactions between warp and weft direction, leading to a shift of the first part of picture frame experimental curves [40,41]. A small equibiaxial strain (0.15%) has been used in the simulation to take into account the influence of the initial tension of the fabric. The yarn width has been used to compare the evolution of the geometry of the unit cell in experimental and simulation tests. The evolution of yarn broadness is shown on Fig. 11b: the width of the yarn obtained by simulation matches the one obtained experimentally at the end of the test, but its decrease does not have exactly the same shape. This behaviour may come from the very low stiffness of compaction behaviour given by the model of Eq. (25). A model with a stiffer initial behaviour could solve this issue. 7. Conclusions A hyperelastic constitutive law has been proposed to describe the mechanical behaviour of the yarns of textile composite reinforcements when they are submitted to large deformations. Physically based invariants have been determined that highlight the different deformation modes identified in the yarn: elongation, compaction, distortion and longitudinal shear. A hyperelastic law has been developed based on these invariants: each deformation mode has been associated a strain energy function which allows describing the internal efforts in the yarn. This approach makes the physical meaning of the model appear clearly. Because some experimental tests on the yarn extracted from the fabric were very difficult, a methodology has been proposed to identify the parameters of this constitutive law from tests at
the scale of the whole reinforcement: uniaxial tension of the yarn and the fabric and biaxial tension of the fabric have been used. An inverse method, based on a Levenberg–Marquardt algorithm, has been used to identify the parameters of compaction and distortion modes. The former constitutive law shows good agreement with experimental biaxial tension test and picture frame test, despite the hypothesis of uncoupling of the different deformation modes. Nevertheless, this assumption may need to be re-examined in other cases. The knowledge of the shear behaviour of the reinforcement has much importance for the simulation of the reinforcement at the scale of the component: it has notably been shown that this behaviour has a major influence on the apparition of wrinkles during the forming process [42]. The obtained possibility of describing the shape of the yarns in the unit cell under deformation is also of great interest in the optimisation of LCM processes: it allows both the knowledge of the mechanical properties of the final component, and the simulation of the injection of the resin inside the reinforcement. Acknowledgements The work reported here has been partly carried out in the scope of the projects MECAFIBRES and LCM3M of the French research agency (ANR). References [1] Advani SG. Flow and rheology in polymeric composites manufacturing. Amsterdam: Elsevier; 1994. [2] Rudd CD, Long AC. Liquid molding technologies. Woodhead Publishing Limited; 1997. [3] Parnas RS. Liquid composite molding. Hanser, Garner Publications; 2000.
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631 [4] Boisse P, Gasser A, Hagège B, Billoet JL. Analysis of the mechanical behaviour of woven fibrous material using virtual tests at the unit cell level. J Mater Sci 2005;40:5955–62. [5] Badel P, Vidal-Sallé E, Maire E, Boisse P. Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale. Compos Sci Technol 2008;68:2433–40. [6] Creech G, Pickett AK. Meso-modelling of non-crimp fabric composites for coupled drape and failure analysis. J Mater Sci 2006;41:6725–36. [7] Zhou G, Sun X, Wang Y. Multi-chain digital element analysis in textile mechanics. Compos Sci Technol 2004;64:239–44. [8] Durville D. A finite element approach of the behaviour of woven materials at microscopic scale. In: 11th Euromech-Mecamat conference. Mechanics of microstructured solids: cellular materials, fiber reinforced solids and soft tissues, Torino, Italy; 2008. [9] Loix F, Badel P, Orgéas L, Geindreau C, Boisse P. Woven fabric permeability: from textile deformation to fluid flow mesoscale simulations. Compos Sci Technol 2008;68:1624–30. [10] 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. Compos Part A: Appl Sci Manuf 1997;28:801–16. [11] Peng X, Cao J. A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics. Compos Part A: Appl Sci Manuf 2005;36:859–74. [12] Willems A. Forming of Textile Composites, PhD Thesis, K.U. Leuven, Belgium; 2008. [13] Aimène Y, Vidal-Sallé E, Hagège B, Sidoroff F, Boisse P. A hyperelastic approach for composite reinforcement large deformation analysis. J Compos Mater 2010;44:5–26. [14] Rivlin RS. Large elastic deformations of isotropic materials. Philos Trans Roy Soc A: Math, Phys Eng Sci 1948;241:379–97. [15] Ogden RW. Non-linear elastic deformations. Paris: John Wiley; 1984. [16] Diani J, Brieu M, Vacherand J-M, Rezgui A. Directional model for isotropic and anisotropic hyperelastic rubber–like materials. Mech Mater 2004;36:313–21. [17] Guo Z, Peng X, Moran B. Large deformation response of a hyperelastic fibre reinforced composite: theoretical model and numerical validation. Compos Part A: Appl Sci Manuf 2007;38:1842–51. [18] Agoras M, Lopez-Pamies O, PonteCastaneda P. A general hyperelastic model for incompressible fiber-reinforced elastomers. J Mech Phys Solids 2009;57:268–86. [19] Holzapfel GA, Gasser TC. A New constitutive framework for arterial wall mechanics and a comparative study of material models. J Elasticity 2000;61:1–48. [20] Hirokawa S, Tsuruno R. Three-dimensional deformation and stress distribution in an analytical/computational model of the anterior cruciate ligament. J Biomech 2000;33:1069–77. [21] Boehler JP. A simple derivation of representations for non-polynomial constitutive equations in some case of anisotropy. J Appl Math Mech 1979;59:157–67. [22] Zhang JM, Rychlewski J. Structural tensors for anisotropic solids. Arch Mech 1990;42:267–77. [23] Zheng QS. Theory of representations for tensor functions – a unified invariant approach to constitutive equations. Appl Mech Rev 1994;47(11):545–86.
1631
[24] Latil P, Orgéas L, Geindreau C, Dumont PJJ, Rolland du Roscoat S. Towards the 3D in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres. Compos Sci Technol 2011;71:480–8. [25] Koissin V, Ivanov DS, Lomov SV, Verpoest I. Fibre distribution inside yarns of textile composite: geometrical and FE modelling. In: Proceedings of the 8th international conference on textile composites (TexComp-8), Nottingham; 2006 [CD edition]. [26] Ciarlet P. Mechanical elasticity: three dimensional elasticity. Elsevier; 1988. [27] Criscione JC, Douglas AS, Hunter WC. Physically based strain invariant set for materials exhibiting transversely isotropic behaviour. J Mech Phys Solids 2001;49:871–97. [28] DeBotton G, Hariton I, Socolsky E. Neo-Hookean fiber-reinforced composites in finite elasticity. J Mech Phys Solids 2006;54:533–59. [29] Buet-Gautier K, Boisse P. Experimental analysis and modelling of biaxial mechanical behaviour of woven composite reinforcements. Exp Mech 2001;41:260–9. [30] Hagège B. Simulation du comportement mécanique des renforts fibreux en grandes transformations. Application aux renforts tricotés, PhD Thesis, Ecole Nationale Supérieure d’Arts et Métiers, France; 2004. [31] Potluri P, Sagar TV. Compaction modelling of textile preforms for composite structures. Compos Struct 2008;86:177–85. [32] Cai Z, Gutowski T. The 3-D deformation behavior of a lubricated fiber bundle. J Compos Mater 1992;26:1207–37. [33] Hivet G, Boisse P. Consistent mesoscopic mechanical behaviour model for woven composite reinforcements in biaxial tension. Compos Part B: Eng 2008;39:345–61. [34] Gorczyca J, Sherwood J, Chen J. A friction model for thermostamping commingled glass–polypropylene woven fabrics. Compos Part A: Appl Sci Manuf 2007;38:393–406. [35] Marquardt DW. An algorithm for least squares estimation of nonlinear parameters. J Soc Ind Appl Math 1963;11:431–41. [36] Schnur DS, Zabaras N. An inverse method for determining elastic material properties and a material interface. Int J Numer Meth Eng 1992;33: 2039–57. [37] Peng XQ, Cao J, Chen J, Xue P, Lussier DS, Liu L. Experimental and numerical analysis on normalization of picture frame tests for composite materials. Compos Sci Technol 2004;64:11–21. [38] Harrison P, Clifford MJ, Long AC. Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments. Compos Sci Technol 2004;64:1453–65. [39] Cao J, Akkerman R, Boisse P, Chen J, et al. Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results. Compos Part A: Appl Sci Manuf 2008;39:1037–53. [40] Launay J, Hivet G, Duong AV, Boisse P. Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements. Compos Sci Technol 2008;68:506–15. [41] Lomov SV, Verpoest I. Model of shear of woven fabric and parametric description of shear resistance of glass woven reinforcements. Compos Sci Technol 2006;66:919–33. [42] Boisse P, Hamila N, Vidal-Sallé E, Dumont F. Simulation of wrinkling during textile composite reinforcement forming, Influence of tensile, in-plane shear and bending stiffnesses. Compos Sci Technol 2011;71:683–92.
Contents lists available at ScienceDirect
Composites Science and Technology journal homepage: www.elsevier.com/locate/compscitech
Hyperelastic modelling for mesoscopic analyses of composite reinforcements A. Charmetant, E. Vidal-Sallé, P. Boisse ⇑ Université de Lyon, LaMCoS, INSA-Lyon, F-69621, France
a r t i c l e
i n f o
Article history: Received 21 April 2011 Received in revised form 24 June 2011 Accepted 3 July 2011 Available online 20 July 2011 Keywords: A. Fabrics/textiles B. Non-linear behaviour C. Anisotropy C. Finite element analysis Hyperelasticity
a b s t r a c t A hyperelastic constitutive law is proposed to describe the mechanical behaviour of fibre bundles of woven composite reinforcements. The objective of this model is to compute the 3D geometry of the deformed woven unit cell. This geometry is important for permeability calculations and for the mechanical behaviour of the composite into service. The finite element models of a woven unit cell can also be used as virtual mechanical tests. The highlight of four deformation modes of the fibre bundle leads to definition of a strain energy potential from four specific invariants. The parameters of the hyperelastic constitutive law are identified in the case of a glass plain weave reinforcement thanks to uniaxial and equibiaxial tensile tests on the fibre bundle and on the whole reinforcement. This constitutive law is then validated in comparison to biaxial tension and in-plane shear tests. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Their good mass-mechanical properties ratio and their good corrosion and fatigue properties make textile composites essential materials when saving mass is a main issue. These materials can be manufactured by Liquid Composite Molding (LCM) processes [1– 3]: the textile reinforcement is preformed in a first stage, then the resin is injected in the preform. The reinforcement is made of continuous fibres grouped into yarns (3000–48,000 fibres in a carbon yarn, 1000–12,000 fibres in a glass yarn). These yarns are combined following textile architectures: woven, knitted braided or non-crimped fabrics. These assemblies give the textile reinforcements a multi-scale nature (Fig. 1). This paper focuses on woven fabrics, but the model which will be proposed for the yarn would be similar for other architectures. The macroscopic scale refers to the whole component level, with dimensions ranging from ten centimetres to several metres. At this level, a woven fabric can be seen as a continuous material whose behaviour is very specific and includes high anisotropy. At the mesoscopic scale, the woven reinforcement is seen as interlaced tows, respectively the warp and the weft (or fill) yarns. Consequently, the working scale corresponds to the yarn dimension, typically one to several millimetres. For periodic materials, mesoscopic models consider the smallest elementary pattern, which can represent the whole fabric by translations. That domain is called the representative unit cell (RUC) [4–6]. Finally, the rein-
⇑ Corresponding author. Tel.: +33 4 72 43 63 96; fax: +33 4 72 43 85 25. E-mail address: [email protected] (P. Boisse). 0266-3538/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compscitech.2011.07.004
forcement can be analysed at the microscopic scale. In this case each fibre is considered as a 3D beam interacting with its numerous neighbours [7,8]. At the microscopic level, the characteristic dimension is about one to several micrometres. This is the only scale at which the material is actually continuous. During the forming process, the woven unit cell reaches its position and geometry in the final composite. This geometry is important regarding the resin flow [9,10] and the final mechanical properties of the composite part. The changes of geometry of the fibre bundles are analysed using mesoscopic scale modelling. The studied domain is then the representative unit cell, which gives the solid skeleton for permeability evaluation. In addition to permeability calculation, mesoscopic analyses can be used as virtual mechanical tests: they allow the mechanical properties of a woven reinforcement to be obtained without performing the experiments and hence without necessarily manufacturing the reinforcement under consideration [4]. In mesoscopic computations the yarn is considered as a continuum but with a very specific behaviour which accounts for his fibrous and discontinuous composition. Some mesoscopic analyses of the woven unit cell submitted to tension, in-plane shear and compression have been set up based on a hypoelastic behaviour of the yarn [4,5,11,12]. These hypoelastic models require the use of an objective derivative, i.e. of a derivative in a frame as close as possible of the fibre orientation. Despite being fairly frequently employed, these models have drawbacks, especially they do not allow a complete recovering after a closed loop loading path. The objective of the development of a hyperelastic constitutive model is to avoid these drawbacks. In addition, a hyperelastic model has been developed recently to analyse the mechanical behaviour of macroscopic textile composite reinforcements modelled as membranes with
1624
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Macroscopic scale Mesoscopic scale Microscopic scale Fig. 1. Different scales for a woven composite reinforcement.
satisfactory results [13]. Nevertheless the model necessary to analyse the deformation of the fibre bundle is quite different since it must be 3D with a single fibre direction while the previous macroscopic model is 2D with two fibre directions. Furthermore the transverse deformation of the yarn (i.e. its compression and shape change), which is essential to the yarn behaviour, cannot be described by a macroscopic membrane approach. An important research effort has been devoted to the development of hyperelastic models through the last decades. Isotropic models have been developed first, to model the behaviour of materials which can be subjected to very large elastic deformation, such as rubber–like materials. These models are usually sorted into two classes, depending on whether the strain energy depends on invariants of the strain tensor in use (e.g. Mooney–Rivlin models [14]), or on principal elongations (e.g. Ogden models [15]). The development of hyperelastic constitutive models for anisotropic materials appeared more recently. They are used to model the mechanical behaviour of strongly anisotropic materials which can undergo large elastic deformations, e.g. reinforced elastomers [16–18] and some biomaterials (biological tissues [19], tendons and ligaments [20]). The formulation of anisotropic hyperelastic models is usually based on structural tensors [21–23] which describe the local orientation (i.e. the preferred directions) of the material. In this paper an anisotropic hyperelastic model is proposed and validated for the mechanical behaviour of textile composite reinforcements. It is based on specific physically motivated invariants. First the specificities of the yarn behaviour are analysed in order to formulate the different assumptions of this work. The mechanical behaviour model of the yarn is then presented, and the identification of its parameters is performed in the case of a glass plain weave, by use of uniaxial and biaxial tensile tests. Finally the model is used to simulate the in-plane shear behaviour of the woven unit cell, and the obtained results are compared to experimental data. 2. Mechanical behaviour of the fibre bundle and assumptions 2.1. Transverse isotropy The yarn is an assembly of fibres oriented approximately in the same direction (Fig. 2). The fibres are assumed to be numerous and compact enough in order to avoid independent movements. In this case, it is possible to consider the yarn as a continuum, and consequently to define a continuous material equivalent to the fibre bundle. In a recent study [24], an X-ray tomography technique was combined to a mechanical test on a fibre bundle. This ap-
Fig. 2. Transverse isotropy of the yarn in a textile reinforcement. X-ray tomography imaging.
proach allowed an accurate observation of the motion of fibres inside the yarn, and showed that such a continuous approach is wellfounded. The homogenized material has a preferred direction which is the fibre one. The spatial distribution of fibres inside the yarn has been analysed on deformed cross sections and, even if this observation may differ for other types of fabrics [25], it has been concluded that this distribution is isotropic [5] for the considered fabric. The ‘‘yarn’’ material will then be assumed to be transversely isotropic.
2.2. Deformation modes A phenomenological approach based on experimental observations is used to define the constitutive law. It is consequently necessary to describe the specific deformation modes of the yarn, which are related to the behaviour of the fibres and to their interactions at the microscopic scale. Four deformation modes are considered: the elongation in the fibre direction, the compaction in the transverse section of the yarn, the distortion (shear) in the transverse section and the shear along the fibre direction (Fig. 3). The last one (longitudinal shear) is a measure of the angle between the covariant cross section and the initial cross section. It is then compound of both longitudinal shear deformations. This deformation mode mainly controls the bending rigidity of the yarn. The experimental analysis of the tensile behaviour of the yarn is quite easy, so its identification can be performed directly from a tensile test on the yarn. This is not the case for the three other deformation modes. The weaving gives cohesion to the fibre bundle in the yarn and it is not easy to define elementary tests for the compaction and shear on a single yarn because of the lack of this cohesion if the yarn is extracted from the fabric. An inverse approach is used in the present work for the determination of the material parameters in transverse compression and shear (modes
1625
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
(a)
(b)
(c)
(d)
Fig. 3. Deformation modes of the yarn: (a) elongation, (b) compaction of the cross section, (c) distortion of the cross section and (d) longitudinal shear. The arrow indicates the direction of fibres.
(b–d) Fig. 3). It is performed on uniaxial and equibiaxial tensile tests because these three modes are present in these cases. 2.3. Uncoupling assumptions The proposed model does not take into account the probably existing couplings between the behaviours corresponding to the deformation modes defined above. It is then assumed that, in the range of deformation simulated and for the considered loading cases, the strain energy function can be split into four parts, each being a function of only one deformation mode. 3. Hyperelastic constitutive equation 3.1. Finite displacement kinematics
~ CQ T Þ; ~ ¼ wðQ wðCÞ
8Q 2 G
ð6Þ
The symmetry group of an isotropic material includes all the rotations and reflections (i.e. G ¼ O3 , the group of isometries of R3 ). Con~ is mathematically isotropic. For an sequently, the strain energy w ~ is mathematically anisotropic. In anisotropic material, GˆO3 so w such a case, it has been shown [21–23] that an isotropic function ^ exists which depends on the right Cauchy–Green strain tensor w C and on N structural tensors Gi defining the symmetry group, such as:
^ ~ ¼ wðC; G1 ; :::; GN Þ wðCÞ
Let X be the initial position of a material point and x its position in the current configuration. v is the transformation from the initial to the current configuration: x ¼ vðXÞ. The deformation gradient tensor is then:
ð7Þ
The representation theorems for transversely isotropic functions of scalars, vectors and/or symmetric tensor variables define a finite ~ such as: number of scalar functions ðI1 ; . . . ; In Þ and a function w _
^ wðC; G1 ; . . . ; GN Þ ¼ wðI1 ; . . . ; In Þ
@x F¼ @X
ð1Þ
This tensor allows the transformation of any vector A in the initial configuration into a vector a of the current (deformed) configuration. The Cauchy–Green deformation tensor C ¼ F T F, from which the variation of the length of A can be deduced, is defined as follows:
kak2 ¼ ðF AÞ ðF AÞ ¼ A ðF T FÞ A ¼ A C A
ð2Þ
The invariants of tensor C are usually defined by:
I1 ¼ traceðCÞ;
1 I2 ¼ ðtraceðCÞ2 traceðC 2 ÞÞ; 2
I3 ¼ detðCÞ
ð3Þ
3.2. Hyperelastic equations
ð4Þ
where F is the deformation gradient tensor and S is the second Piola–Kirchhoff stress tensor. To ensure that this constitutive equation fulfils the principle of ~ of the right Cauchy Green material frame indifference, a function w strain tensor C must exist as [26]:
~ ~ T FÞ ¼ wðCÞ ¼ wðFÞ wðF
Number n depends on the number and type of arguments of the isotropic function under consideration. Finally, the strain energy function will be defined in the form (8). For sake of simplicity, it will be denoted w in the sequel. The symmetry group of a transversely isotropic material is characterised, in the initial configuration, by an unit vector M which represents the preferred direction. This allows the definition of a structural tensor M ¼ M M. Then, the strain energy potential becomes:
ð9Þ
where I1, I2, I3 are the standard invariants of C defined in (3). II4 and II5 are mixed invariants defined from the structural tensor M by:
II4 ¼ C : M
and II5 ¼ C 2 : M
ð10Þ
The second Piola–Kirchhoff stress tensor is then defined as:
In order to define a hyperelastic material, an elastic strain energy potential per unit volume w must exist which only depends on the current strain state [26]. As the behaviour is reversible, it does not dissipate energy and naturally fulfils the Clausius–Duhem inequality. These assumptions allow the constitutive equation of a hyperelastic material to be written as:
@wðFÞ @C
ð8Þ
w ¼ wðI1 ; I2 ; I3 ; II4 ; II5 Þ
The local change of volume during deformation is described by the pffiffiffiffi Jacobian J ¼ I3 ¼ detðFÞ.
SðFÞ ¼ 2
Consequently, the strain energy potential can be defined directly, ~ using wðCÞ. Moreover, it must satisfy the principle of material symmetry, i.e. the strain energy must be invariant under any transformation belonging to the symmetry group G of the material:
ð5Þ
S¼2
@w @C
@w @I1 @w @I2 @w @I3 @w @II4 @w @II5 þ þ þ þ ¼2 @I1 @C @I2 @C @I3 @C @II4 @C @II5 @C
! ð11Þ
4. Constitutive equation for a yarn In this section, a hyperelastic model is proposed which describes the behaviour of the yarn in textile composite reinforcements. First, ‘‘physically based’’ invariants are defined which allow the description of the different deformation modes of the yarn introduced in Section 2. Then, each mode will be associated a strain energy function from which the stresses in the material can be derived.
1626
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
4.1. Invariants and their corresponding deformation modes An approach based on a decomposition of the deformation gradient [27] is used to determine physically based invariants of the transformation for each deformation mode of the fibre bundle. An orthonormal basis B ¼ fM; N 1 ; N 2 g is defined which allows writing the components of the deformation gradient F of transversely isotropic materials into the form:
2
fm 6 ½FB ¼ 4 0 0
fm1
fm2
3
0
0
0
32
1
½F elong B
2
½F comp B
1
0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 f11 =f22 40
32
0
1
fm1 fm
fm2 fm
3
76 7 0 5 40 1 0 5 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 0 f22 =f11 0 0 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} ½F dist B
ð13Þ
½F sh B
where F elong , F comp , F dist and F sh respectively describe the elongation of the yarn along the fibres direction, the compaction and distortion of the yarn in the transverse direction, and the longitudinal shear of the yarn. Each part of the decomposition can be characterized by a single quantity, except F sh which requires two quantities. Finally, five parameters are needed to describe the deformation with this approach:
aelong ¼ fm ; acomp ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f2 þ f2 ash ¼ m1 2 m2 ; fm
pffiffiffiffiffiffiffiffiffiffiffi f11 f22 ;
tanðcÞ ¼
adist ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi f11 =f22 ;
fm1 fm2
ð14Þ
This is consistent with the fact that five invariants are needed to write the strain energy potential (9) of hyperelastic transversely isotropic materials. With these notations, the components of the transformation gradient can be written as:
2
sffiffiffiffiffiffiffiffiffiffiffiffiffi I5 1 I24
ð17Þ
These invariants are used in the sequel to define the strain energy functions associated to each deformation mode:
ð18Þ
4.2. Strain energy functions associated to each deformation mode
3 0 0 p ffiffiffiffiffiffiffiffiffiffiffi 6 76 0 7 f11 f22 ½FB ¼ 4 0 1 0 5 4 0 5 pffiffiffiffiffiffiffiffiffiffiffi 0 0 1 0 0 f11 f22 |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} fm
Ish ¼
ð12Þ
As previously mentioned, M is a unit vector which direction is the preferred direction of the material (i.e. the direction of fibres) in the initial configuration. A multiplicative split of the matrix (12) can be made in order to separate the parts corresponding to each deformation mode of the yarn:
2
Idist
wðCÞ ¼ wðIelong ; Icomp ; Idist ; Ish Þ
7 0 5 f22
f11
1 1 I3 ; lnðI4 Þ; Icomp ¼ ln 2 4 I4 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 2 1 I 1 I 4 I5 I1 I 4 I 5 pffiffiffiffiffiffiffiffi ¼ ln @ pffiffiffiffiffiffiffiffi þ 1A; 2 2 I3 I4 2 I3 I4
Ielong ¼
3
aelong ash aelong cosðcÞ ash aelong sinðcÞ 7 acomp adist 0 5 0 acomp =adist
6 ½FB ¼ 4 0 0
sffiffiffiffiffiffiffiffiffiffiffiffiffi I ash ¼ 52 1 I4
nl w elong ðIelong Þ if welong ðIelong Þ ¼ wlin ðIelong Þ if elong
Ielong 6 I0elong
ð19Þ
Ielong > I0elong
The strain energy functions corresponding to the non-linear and linear parts of the elongation behaviour of the yarn are respectively defined by the following equations: 2 3 wnl elong ðIelong Þ ¼ a0 þ a1 I elong þ a2 Ielong þ a3 Ielong
ð20Þ
2 wlin elong ðIelong Þ ¼ b0 þ b1 I elong þ b2 I elong
ð21Þ
Introducing the initial yarn section S0 (assumed to be constant along the yarn), its initial stiffnesses K 0elong , its stiffness in the linear part Kelong, and the continuity condition between the two functions lin wnl elong et welong , leads to:
wnl elong ¼
wlin elong ¼
ð15Þ
K 0elong 2 K elong K 0elong 3 Ielong þ Ielong 2S0 6S0 I0elong
ð22Þ
K elong K 0elong 0 2 K elong K 0elong 0 Ielong Ielong Ielong 6S0 2S0 K elong 2 þ I 2S0 elong
ð23Þ
Using (11), the second Piola–Kirchhoff stress tensor is given by:
The quantity c is a measure of the coupling between the two shear modes, namely distortion and longitudinal shear [27,28]. Under the uncoupling assumptions formulated before, the behaviour of the yarn is then assumed not to depend on this quantity. The standard invariants (9) can then be expressed as functions of the quantities (14). Inverting the obtained expressions leads to:
vffiffiffiffiffiffiffiffi ffi usffiffiffiffi u I pffiffiffiffi aelong ¼ I4 ; acomp ¼ t 3 ; I4 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 uI I I I1 I4 I5 t14 5 p pffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffi ffi þ adist ¼ 1; 2 I3 I4 2 I3 I4
4.2.1. Elongation strain energy A tensile test on a single yarn allows the identification of its tensile behaviour. Such a test highlights a response of the yarn in tension which can be close to a linear behaviour (see Fig. 4). Nevertheless, the slight misalignment of the fibres inside the tow can induce an initial non linearity: as the fibres do not have exactly the same length, the overall tension of the yarn corresponds to a non simultaneous tension of the fibres, and consequently to a progressively stiffening behaviour of the yarn [29]. In order to describe both linear and non-linear parts of the response of the yarn in tension, the corresponding strain energy potential is defined piecewise:
ð16Þ
Amongst the previous quantities, only ash vanishes when the material is not deformed (i.e. when F ¼ I). After normalizing the other quantities, the following set of invariants is defined:
Selong
0 K elong K 0elong 2 K elong 1 S0 Ielong þ 2S0 I0elong Ielong ¼ M I4 K elong K 0elong 0 K I þ elong Ielong 2S0
elong
S0
if
Ielong 6 I0elong
if
Ielong > I0elong
ð24Þ
Four parameters have to be identified: S0, I0elong , K 0elong and Kelong. The first one will be determined by measuring the geometrical characteristics of the fabric. To avoid stiffness errors, S0 must exactly match the initial cross-sectional area defined in the finite element model. The three remaining parameters will be identified using the result of the tensile test. 4.2.2. Compaction strain energy The strain energy associated to the transverse compaction of the yarn is linked to several physical phenomena at the microscopic scale (mostly bending and rearrangement of fibres) which
1627
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Identified parameters for elongation
200
S0 (mm²) I elong 0 0 (N) K elong K elong (N)
150
Simulation
Effort (N)
1.0496 0.00134 11660 36714
Experience
100 50 0 0
0.002
0.004
0.006
Deformation of the yarn Fig. 4. Identification of yarn elongation behaviour.
are difficult to study separately. As in [4,30–32], a power-based strain energy function is then proposed:
p K comp jIcomp j wcomp ðIcomp Þ ¼ 0
if
Icomp 6 0
if
Icomp > 0
ð25Þ
ð26Þ
The parameters to be identified are then Kcomp and p. As a direct experimental characterization of this behaviour is difficult, an inverse method is set up, based on equibiaxial tensile tests, to identify both parameters. An equibiaxial tensile test is used rather than a crushing test because the later has been found to be very sensitive to the frictional coefficient between the crushing tools and the fabric, which may make the identification less reliable. Moreover, crushing simulations computation time is much larger, which significantly increases the duration of the identification procedure. 4.2.3. Distortion strain energy As well as those of the compaction behaviour, the parameters of the behaviour of the distortion in the transverse section can not be determined by a direct experiment. In the present work, it has been assumed that the distortion stiffness is constant which leads to the following strain energy function:
wdist ðIdist Þ ¼
1 K dist I2dist 2
ð27Þ
Constant Kdist will be identified at the same time as the parameters of Eq. (25), using the same inverse method. The second Piola–Kirchhoff stress tensor associated to this distortion strain energy is then given by:
Sdist ¼ 2K dist Idist
2I4 I ðI1 I4 I5 ÞC 1 þ I1 þ II54 M 2ðC M þ M CÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 ðI1 I4 I5 Þ2 4I3 I4
1 K sh I2sh 2
ð29Þ
The second Piola–Kirchhoff stress tensor is then given by:
It is assumed that no energy is required for yarn section expansion. It is also assumed that during the forming of a textile composite reinforcement, the yarn compaction stiffness never reaches the bulk stiffness of glass (i.e. the voids in the yarn do not collapse altogether). The second Piola–Kirchhoff stress tensor, corresponding to the yarn compaction, derives from Eq. (25):
p 1 Scomp ðIcomp 6 0Þ ¼ K comp jIcomp jp1 C 1 M 2 I4
wsh ðIsh Þ ¼
ð28Þ
4.2.4. Longitudinal shear strain energy The longitudinal shear strain energy of the yarn is linked to the sliding of its constitutive fibres in their preferred direction. A linear elastic behaviour will be assumed for this deformation mode, which induces the following strain energy potential:
" # 1 1 2I5 Ssh ¼ K sh 2 ðC M þ M CÞ 3 M 2 I4 I4
ð30Þ
During a tensile test on a textile composite reinforcement, this deformation mode is activated: during the non-linear part of the tensile curve that corresponds to a change of undulation of the yarns, before reaching the maximum tensile stiffness, longitudinal shear and elongation are the only working modes, whereas in the linear part, only the longitudinal elongation works (Fig. 7). The identification of Ksh is then possible using an inverse method on the uniaxial tensile test as soon as the longitudinal elongation behaviour is known. 4.3. Combination of all contributions Finally, the whole constitutive equation is the summation of all contributions discussed before:
S¼2
@welong @Ielong @wcomp @Icomp @wdist @Idist @wsh @Ish þ þ þ @Ielong @C @Icomp @C @Idist @C @Ish @C
!
¼ Selong þ Scomp þ Sdist þ Ssh ð31Þ Table 1 summarizes the parameters which have to be identified for each deformation mode, as well as the experimental test used for their identification. The proposed constitutive equation has been implemented as a user material subroutine VUMAT in the ABAQUS/Explicit finite element code. 4.4. Studied reinforcement and finite element model of the yarn For sake of simplicity, the parameter identification has been realized on a glass plain weave which main geometrical parameters are summarized in Table 2. The plain weave presents the great advantage to have a rather simple unit cell: it is possible to obtain a realistic geometrical model of the yarn, using only straight lines and circular arches, avoiding yarn interpenetrations [33]. Fig. 5 shows the geometry used in the present work. The elementary unit cell has several symmetry axes and planes. When the boundary conditions imposed to the unit cell are compatible with those symmetries, it becomes possible to reduce the pattern to be studied up to 1/8 of the unit cell (Fig. 6). For example, it is the case of biaxial tension and crushing. Nevertheless the whole unit cell has to be used in order to simulate in-plane shear. Contact and frictional sliding between the tows is treated using Coulomb friction with a coefficient of 0.3 identified in [34].
1628
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Table 1 Summary of the parameters and their identification approaches. Deformation mode
Elongation
Compaction
Number of parameters Identification mode
4 Geometric observations, tensile test on the yarn
2 1 Equibiaxial tensile test, iterative approach
Plain weave
Yarn width (mm)
Warp: 3.2 Weft: 3.1
Yarn density (yarn/mm)
Warp: 0.251 Weft: 0.248
Crimp (%)
Warp: 0.5 Weft: 0.54
Surface weight (g/m2)
600
-3
Weave type
Energies (mJ.mm )
Table 2 Geometrical characteristics of the studied glass plain weave.
Distortion
Longitudinal shear 1 Tensile test on the fabric
1.00E+00 1.00E-01 1.00E-02 1.00E-03 1.00E-04 1.00E-05 1.00E-06 1.00E-07 1.00E-08 1.00E-09 1.00E-10 1.00E-11 1.00E-12
Welong Wsh Wcomp Wdist
0
0.02
0.04
0.06
0.08
0.1
0.12
Deformation of the fabric Fig. 7. Evolution of strain energies associated to each deformation mode during an uniaxial tensile test on the reinforcement.
(a)
(b)
Circle arches
Straight segments
Circle arches
Fig. 5. Real geometry obtained using X-ray tomography imaging: (a) and corresponding model and (b) for a unit cell of plain weave.
yarn bending). Once the parameters of the elongation behaviour of the yarn have been identified, it becomes, then, possible to identify the parameter Ksh of Eqs. (29) and (30). This method allows determining the value Ksh = 3 MPa. On Fig. 7, the strain energies associated to each deformation mode during a uniaxial tensile test on the fabric are plotted. Since the energies of elongation and longitudinal shear are much greater than those associated to compaction and distortion, this method of identification using this kind of test seems well-founded.
5. Parameters identification
5.3. Identification of compaction and distortion behaviours
5.1. Identification of the behaviour in elongation
As mentioned before, the parameters Kcomp, Kdist and p from Eqs. (25) and (27) are then identified using an inverse method: they are optimised so that the simulated behaviour of the unit cell matches its experimentally observed behaviour. To ensure an accurate identification, the experimental test of reference must mainly depend on these deformation modes of the yarn. In this paper the experimental test of reference chosen is equibiaxial tension which, as seen afterwards, satisfies the former condition. An equibiaxial tension test consists in imposing the same elongation to both warp and weft directions of the woven reinforcement. Such a test induces a significant amount of compaction and distortion: because the tensile stiffness of the yarn is very high, a fabric which undergoes equibiaxial tension starts adapting itself to the deformation by reducing its crimp (i.e. the undulations of the yarns). This crimp reduction is made possible by both compac-
The initial section of the yarn, S0, is given by mesoscale observations of the fabric, using X-ray tomography for instance. The three other parameters, namely I0elong , K 0elong and Kelong, are identified thanks to a tensile test on a single yarn. As the yarn may have been damaged during the weaving, this test must be performed on a yarn extracted from the fabric rather than a yarn before weaving. The parameters obtained by this method for the glass fibre yarn, as well as the experimental and identified curves, are shown on Fig. 4. 5.2. Identification of the longitudinal shear behaviour During an uniaxial tensile test on the fabric, the main deformation modes of the yarn are elongation and longitudinal shear (i.e.
(a)
(b)
(c)
Fig. 6. (a) Plain weave. (b) Elementary unit cell for in-plane shear. (c) 1/8 of the unit cell for biaxial tension.
1629
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
Identified parameters for compaction
140
K comp (N) pcomp
120
Experience (warp direction)
Identified parameter for distortion K dist 0.06
Effort per yarn (N)
17800 7.72
Experience (weft direction) Initial simulation (10000,5,0.05)
100
Loop 1 (9300,6.82,0.0542) Loop 2 (10300,7.16,0.0561)
80
Loop 3 (17800,7.72,0.0608)
60 40 20 0 0
0.001 0.002 0.003 0.004 0.005 0.006 Deformation of the fabric
Fig. 8. Optimisation of the parameters of compaction and distortion using equibiaxial tension.
200 180
Effort per yarn (N)
tion and distortion. That is why the non-linear part of the curve of equibiaxial tension (Fig. 8) is representative of the behaviour of the yarn in its transverse plane. The inverse method which has been set up consists in minimizing a least square estimation of the error between the experimental curve of equibiaxial tension, and the curve obtained by simulation. An iterative Levenberg–Marquardt algorithm [35,36] has been used to determine the three parameters of compaction and distortion. Fig. 8 illustrates the successive curves obtained during this optimisation. The proposed method allows to obtain a result of simulation very close to the observed behaviour, with only few iterations. Nevertheless, the existence of local minima of the error estimation induces that the initial parameters must be close enough to the final solution so that the algorithm converges. This imposes a first trial and error approach before starting the optimisation procedure.
Experience Simulation
160
Yarn k=0.5 free
140 120
k=1
100 80
k=2
60 40 20 0 0
0.002
0.004
0.006
0.008
Deformation of the fabric Fig. 9. Biaxial tension curves obtained after identification of the parameters.
6. Validation of the hyperelastic law In this section, the previous hyperelastic constitutive law and the identification of its parameters are validated: the parameters previously identified are used to perform simulations in other cases and the results are compared to experimental values. A first validation is established by comparing simulation and experimental results of biaxial tension with different warp/weft elongation ratios. Fig. 9 presents the whole network of experimentally and numerically obtained curves. Remind that the curves ‘‘yarn’’, ‘‘free’’ and ‘‘k = 1’’ have been used for identification. The validation can then be performed by comparing simulation and experimental curves for the cases ‘‘k = 0.5’’ and ‘‘k = 2’’. The two numerically obtained curves show good agreement as regard as experimental curves, which allows a first validation of the proposed hyperelastic constitutive law. Biaxial tension tests with k = 0.5 and k = 2 elongation ratios are, kinematically, relatively close to the k = 1 test used for the identification. That is why another validation is performed using an experimental test which kinematics is very different from biaxial tension. Then, a simulation of the shear behaviour of the unit cell has been performed and compared to the results of the picture frame test [37–40]. In order to compare the meso-scale simulation to the macro-scale picture frame test, it is assumed that the fabric and the unit cell have the same mechanical behaviour, i.e. that the shear deformation in the fabric during a picture frame test is uniformly distributed. The periodicity and boundary conditions imposed to the unit cell in this simulation of shear are similar to
those described in [5]. The initial and deformed configurations of the unit cell are shown on Fig. 10. Two criteria are used to compare the results: one which allows validating the overall efforts induced in the unit cell and one to validate its deformed geometry. In order to compare the results obtained by different measurement devices, or to compare the inplane shear behaviour of different types of fabrics, a quantity must be defined which characterize the efforts in the unit cell. Different quantities have been introduced in previous papers. In the sequel the torque Cs per unit initial area introduced in [40] is used:
Cs ¼ F
cosðcÞ L2fabric 2cosðp=4 c=2Þ
Lframe
ð32Þ
with F the effort measured on the frame, Lframe the side of the frame, Lfabric the side of the square of fabric inside the frame and c the angle variation between warp and weft directions (i.e. the shear angle). On the other hand, an energy-based approach [40] allows the calculation of this torque per unit initial area (32) from the simulation:
Cs ¼
_ W Su c_
ð33Þ
_ is the time-derivative of the external work, Su the initial where W surface of the unit cell and c_ the time-derivative of the shear angle. The comparison of the torque obtained from experiments and simulations is presented on Fig. 11a. On this figure, the shear torques computed during the simulation and measured during the experi-
1630
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631
0.7 Picture frame test
0.6
Simulation
0.5 0.4 0.3 0.2 0.1 0 0
(a)
20
40
60
Shear angle (deg)
Relative broadness of the yarns
Torque per unit initial area (N)
Fig. 10. Unit cell shear simulation: (a) undeformed configuration and (b) deformed configuration.
1 0.8 0.6 0.4
Experiment
0.2
Simulation
0 0
10
(b)
20
30
40
50
60
Shear angle (deg)
Fig. 11. (a) Evolution of normalized shear force during numerical and experimental tests. (b) Evolution of yarn broadness during numerical and experimental tests.
ment are close. During the preparation of shear samples for experiments, a slight tension of the fabric is required to ensure a good positioning of the fabric in the picture frame. This tension induces stronger interactions between warp and weft direction, leading to a shift of the first part of picture frame experimental curves [40,41]. A small equibiaxial strain (0.15%) has been used in the simulation to take into account the influence of the initial tension of the fabric. The yarn width has been used to compare the evolution of the geometry of the unit cell in experimental and simulation tests. The evolution of yarn broadness is shown on Fig. 11b: the width of the yarn obtained by simulation matches the one obtained experimentally at the end of the test, but its decrease does not have exactly the same shape. This behaviour may come from the very low stiffness of compaction behaviour given by the model of Eq. (25). A model with a stiffer initial behaviour could solve this issue. 7. Conclusions A hyperelastic constitutive law has been proposed to describe the mechanical behaviour of the yarns of textile composite reinforcements when they are submitted to large deformations. Physically based invariants have been determined that highlight the different deformation modes identified in the yarn: elongation, compaction, distortion and longitudinal shear. A hyperelastic law has been developed based on these invariants: each deformation mode has been associated a strain energy function which allows describing the internal efforts in the yarn. This approach makes the physical meaning of the model appear clearly. Because some experimental tests on the yarn extracted from the fabric were very difficult, a methodology has been proposed to identify the parameters of this constitutive law from tests at
the scale of the whole reinforcement: uniaxial tension of the yarn and the fabric and biaxial tension of the fabric have been used. An inverse method, based on a Levenberg–Marquardt algorithm, has been used to identify the parameters of compaction and distortion modes. The former constitutive law shows good agreement with experimental biaxial tension test and picture frame test, despite the hypothesis of uncoupling of the different deformation modes. Nevertheless, this assumption may need to be re-examined in other cases. The knowledge of the shear behaviour of the reinforcement has much importance for the simulation of the reinforcement at the scale of the component: it has notably been shown that this behaviour has a major influence on the apparition of wrinkles during the forming process [42]. The obtained possibility of describing the shape of the yarns in the unit cell under deformation is also of great interest in the optimisation of LCM processes: it allows both the knowledge of the mechanical properties of the final component, and the simulation of the injection of the resin inside the reinforcement. Acknowledgements The work reported here has been partly carried out in the scope of the projects MECAFIBRES and LCM3M of the French research agency (ANR). References [1] Advani SG. Flow and rheology in polymeric composites manufacturing. Amsterdam: Elsevier; 1994. [2] Rudd CD, Long AC. Liquid molding technologies. Woodhead Publishing Limited; 1997. [3] Parnas RS. Liquid composite molding. Hanser, Garner Publications; 2000.
A. Charmetant et al. / Composites Science and Technology 71 (2011) 1623–1631 [4] Boisse P, Gasser A, Hagège B, Billoet JL. Analysis of the mechanical behaviour of woven fibrous material using virtual tests at the unit cell level. J Mater Sci 2005;40:5955–62. [5] Badel P, Vidal-Sallé E, Maire E, Boisse P. Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale. Compos Sci Technol 2008;68:2433–40. [6] Creech G, Pickett AK. Meso-modelling of non-crimp fabric composites for coupled drape and failure analysis. J Mater Sci 2006;41:6725–36. [7] Zhou G, Sun X, Wang Y. Multi-chain digital element analysis in textile mechanics. Compos Sci Technol 2004;64:239–44. [8] Durville D. A finite element approach of the behaviour of woven materials at microscopic scale. In: 11th Euromech-Mecamat conference. Mechanics of microstructured solids: cellular materials, fiber reinforced solids and soft tissues, Torino, Italy; 2008. [9] Loix F, Badel P, Orgéas L, Geindreau C, Boisse P. Woven fabric permeability: from textile deformation to fluid flow mesoscale simulations. Compos Sci Technol 2008;68:1624–30. [10] 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. Compos Part A: Appl Sci Manuf 1997;28:801–16. [11] Peng X, Cao J. A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics. Compos Part A: Appl Sci Manuf 2005;36:859–74. [12] Willems A. Forming of Textile Composites, PhD Thesis, K.U. Leuven, Belgium; 2008. [13] Aimène Y, Vidal-Sallé E, Hagège B, Sidoroff F, Boisse P. A hyperelastic approach for composite reinforcement large deformation analysis. J Compos Mater 2010;44:5–26. [14] Rivlin RS. Large elastic deformations of isotropic materials. Philos Trans Roy Soc A: Math, Phys Eng Sci 1948;241:379–97. [15] Ogden RW. Non-linear elastic deformations. Paris: John Wiley; 1984. [16] Diani J, Brieu M, Vacherand J-M, Rezgui A. Directional model for isotropic and anisotropic hyperelastic rubber–like materials. Mech Mater 2004;36:313–21. [17] Guo Z, Peng X, Moran B. Large deformation response of a hyperelastic fibre reinforced composite: theoretical model and numerical validation. Compos Part A: Appl Sci Manuf 2007;38:1842–51. [18] Agoras M, Lopez-Pamies O, PonteCastaneda P. A general hyperelastic model for incompressible fiber-reinforced elastomers. J Mech Phys Solids 2009;57:268–86. [19] Holzapfel GA, Gasser TC. A New constitutive framework for arterial wall mechanics and a comparative study of material models. J Elasticity 2000;61:1–48. [20] Hirokawa S, Tsuruno R. Three-dimensional deformation and stress distribution in an analytical/computational model of the anterior cruciate ligament. J Biomech 2000;33:1069–77. [21] Boehler JP. A simple derivation of representations for non-polynomial constitutive equations in some case of anisotropy. J Appl Math Mech 1979;59:157–67. [22] Zhang JM, Rychlewski J. Structural tensors for anisotropic solids. Arch Mech 1990;42:267–77. [23] Zheng QS. Theory of representations for tensor functions – a unified invariant approach to constitutive equations. Appl Mech Rev 1994;47(11):545–86.
1631
[24] Latil P, Orgéas L, Geindreau C, Dumont PJJ, Rolland du Roscoat S. Towards the 3D in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres. Compos Sci Technol 2011;71:480–8. [25] Koissin V, Ivanov DS, Lomov SV, Verpoest I. Fibre distribution inside yarns of textile composite: geometrical and FE modelling. In: Proceedings of the 8th international conference on textile composites (TexComp-8), Nottingham; 2006 [CD edition]. [26] Ciarlet P. Mechanical elasticity: three dimensional elasticity. Elsevier; 1988. [27] Criscione JC, Douglas AS, Hunter WC. Physically based strain invariant set for materials exhibiting transversely isotropic behaviour. J Mech Phys Solids 2001;49:871–97. [28] DeBotton G, Hariton I, Socolsky E. Neo-Hookean fiber-reinforced composites in finite elasticity. J Mech Phys Solids 2006;54:533–59. [29] Buet-Gautier K, Boisse P. Experimental analysis and modelling of biaxial mechanical behaviour of woven composite reinforcements. Exp Mech 2001;41:260–9. [30] Hagège B. Simulation du comportement mécanique des renforts fibreux en grandes transformations. Application aux renforts tricotés, PhD Thesis, Ecole Nationale Supérieure d’Arts et Métiers, France; 2004. [31] Potluri P, Sagar TV. Compaction modelling of textile preforms for composite structures. Compos Struct 2008;86:177–85. [32] Cai Z, Gutowski T. The 3-D deformation behavior of a lubricated fiber bundle. J Compos Mater 1992;26:1207–37. [33] Hivet G, Boisse P. Consistent mesoscopic mechanical behaviour model for woven composite reinforcements in biaxial tension. Compos Part B: Eng 2008;39:345–61. [34] Gorczyca J, Sherwood J, Chen J. A friction model for thermostamping commingled glass–polypropylene woven fabrics. Compos Part A: Appl Sci Manuf 2007;38:393–406. [35] Marquardt DW. An algorithm for least squares estimation of nonlinear parameters. J Soc Ind Appl Math 1963;11:431–41. [36] Schnur DS, Zabaras N. An inverse method for determining elastic material properties and a material interface. Int J Numer Meth Eng 1992;33: 2039–57. [37] Peng XQ, Cao J, Chen J, Xue P, Lussier DS, Liu L. Experimental and numerical analysis on normalization of picture frame tests for composite materials. Compos Sci Technol 2004;64:11–21. [38] Harrison P, Clifford MJ, Long AC. Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments. Compos Sci Technol 2004;64:1453–65. [39] Cao J, Akkerman R, Boisse P, Chen J, et al. Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results. Compos Part A: Appl Sci Manuf 2008;39:1037–53. [40] Launay J, Hivet G, Duong AV, Boisse P. Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements. Compos Sci Technol 2008;68:506–15. [41] Lomov SV, Verpoest I. Model of shear of woven fabric and parametric description of shear resistance of glass woven reinforcements. Compos Sci Technol 2006;66:919–33. [42] Boisse P, Hamila N, Vidal-Sallé E, Dumont F. Simulation of wrinkling during textile composite reinforcement forming, Influence of tensile, in-plane shear and bending stiffnesses. Compos Sci Technol 2011;71:683–92.