Le présent travail a pour objectif de présenter une étude de sensibilité des modèles numériques réalisés avec le code ABAQUS, vis-à-vis de la variation des paramètres du procédé d'emboutissage, ainsi que l'effet d'orientation initiale du renfort, le type (coque ou membrane) et la taille du maillage. La simulation de la mise en forme est réalisée à l'échelle macroscopique en considérant le renfort comme un milieu continu. Cette approche continue s'appuie sur une loi de comportement hypoélastique qui a l'aptitude de suivre la rotation des fibres (directions d'anisotropie) au cours de la mise en forme. Cette loi de comportement est implémentée dans le code de calcul des éléments finis ABAQUS /explicit en utilisant une subroutine VUMAT

Arfaoui, Makrem

Aridhi, Abderrahmen

Boisse, Philippe

Mabrouki, Tarek

Naouar, Naim

Zarroug, Malek

I-Revues

23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
 
Modeling and numerical simulation of the drawing 
of a dry woven reinforcement: Sensitivity of the 
shear angle to the parameters of the process 
Modélisation et simulation numérique de l’emboutissage d’un renfort 
tissé sec : Sensibilité de l’angle de cisaillement aux paramètres du 
procédé 
A. ARIDHI a,b, T. MABROUKI a, M. ARFAOUI a, N. NAOUAR b  
P. BOISSE b  M. ZARROUG c  
a Université de Tunis El Manar, École Nationale d’Ingénieurs de Tunis, LR-11-ES19 Laboratoire 
de Mécanique Appliquée et Ingénierie, 1002 Tunis, Tunisie 
b Laboratoire de Mécanique des Contacts et des Solides UMR CNRS 5514, INSA de Lyon, Bâtiment 
Jacquard, 27 Rue Jean Capelle, 69621 Villeurbanne Cedex, France 
c 
PSA Peugeot-Citroën – Direction Scientifique – DRD/DRIA/DSTF/SMMS 
E-mail : abderrahmen.aridhi@insa-lyon.fr 
 
Résumé : (16 gras) 
Le présent travail a pour objectif de présenter une étude de sensibilité des modèles numériques 
réalisés avec le code ABAQUS, vis-à-vis la variation des paramètres du procédé d’emboutissage,  
l'effet d'orientation initiale du renfort, le type (coque ou membrane) et la taille du maillage. La 
simulation de la mise en forme est réalisée à l’échelle macroscopique en considérant le renfort comme 
un milieu continue. Cette approche continue s’appuie sur une loi de comportement hypo-élastique non 
orthogonale qui a l’aptitude à suivre la rotation des fibres (direction d’anisotropie) au cours de la 
mise en forme. Cette loi de comportement est implémentée dans le code de calcul des éléments finis 
ABAQUS/explicit en utilisant une subroutine VUMAT. 
 
Abstract : (16 gras) 
The present work aims to present a sensitivity study of the numerical models that are modelled with 
ABAQUS code, towards the variation of the process parameters, the initial orientation of fabric, type 
(shell or membrane) and size of mesh. Forming simulation has been realized at the macroscopic scale 
by considering the fabric as a continuum medium. The mechanical approach is based on non-
orthogonal hypo-elastic behaviour that can track the rotation of fiber during the forming process. The 
constitutive model has been implemented in a commercial FE code (ABAQUS Explicit) via a user 
material subroutine VUMAT. 
 
Keywords: hypo-elasticity / Continuous Approach / Textile Reinforcement 
 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
1. Introduction 
Continuous fiber reinforced thermoplastic composites have been increasingly employed in several 
industrial fields such as aerospace, automotive, energy areas, defense and sports, due to enhanced and 
multi-functional material properties as compared to sheet metals, for instance their high specific 
strength and stiffness [1].  In recent years, the automakers are faced with two difficulties. The first one 
is the reduction of polluting emissions by decreasing the cars weight and the second challenge is about 
improving the cadence of production and mechanicals properties of the part [2]. In this context, the use 
of composite woven continuous glass fibers and thermoplastic matrix shaped by thermo-stamping, 
constitute the best compromise mechanical performance, cost and cadence of production [2]. Indeed 
the woven composite possess a good formability which results from the stretching of fibers in fiber 
direction and large amount of extension in diagonal direction. According these important 
characteristics, the textile composite can drape into fairly complex shape. In addition, the thermo-
stamping process has a rapid cycle of fabrication (only seconds) regardless of the size and the 
complexity of the parts compared other process (LCM) which are very expensive with long time cycle 
[3]. The aim of this paper is to present a non-orthogonal hypo-elastic model that can track the fiber 
rotation during the forming step then to present a sensitivity study of the numerical models that are 
modelled with ABAQUS code, towards the variation of the process parameters  
2. Computational implementation 
2.1. Hypoelastic formulation: Constitutive model for woven 
fabrics 
The hypo-elastic laws are widely used in FE codes to model the composite's behavior at large strain 
[4][5]: 
 :C D     (1) 
Where
 
and D
 
are, respectively, the Cauchy stress and the strain rate tensors. 

 
is the objective 
derivative of the stress tensor. It can be seen as the time derivative for an observer which is fixed with 
regards to the material. The objective derivative is used to avoid stress change due to rigid body 
rotations [3][5]. It is written as:   
  . .t tdQ Q Q Q
dt
 
 
  
 
  (2) 
In literature, there are several objective derivatives that the main object is avoiding the stress change 
under rigid body rotations. Commercial FE codes, such as ABAQUS use generally two types of 
derivative objective [7]  
- The Green-Naghdi GN stress rate, based on polar rotation, is used in ABAQUS/Explicit 
- The Jauman stress rate, based on corotational frame rotation, is used in ABAQUS/Standard  
 
In order to follow the rotation of fiber yarns in woven fabrics during deformation, we construct two 
vectors which are coincided with the current directions of warp and weft. In reality, the two vectors 
are no orthogonal and no unit but for simplicity they are chosen to be unit vectors. In the current 
configuration, the fiber directions are obtained from the deformation gradient tensor F . By default the 
stress and strain are expressed in the GN base. So it must define the vectors  1 2,e e  that defined the 
GN basis.  In the current configuration,  1 2,e e   are obtained using orthogonal rotation tensor R .  
 
0 0
0 0
. .
;
. .
k k
k k
k k
F e R e
f e
F e R e
    (3) 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
   
 
Figure 1: Schematic of a deformed plain weave structure with shear deformation 
Initially, the initial vectors of the two bases (Fiber and GN) are assumed coincided. Let   denotes the 
angle between the two current vectors of fiber frames and   defines the angle between the two 
independent fiber frames and the orthogonal GN axes (see figure1). To take into account the non-
orthogonally of woven fabrics after forming the dry woven composite, we need to give the expression 
of the stress and deformation tensors in the non-orthogonal frame. In order to express the stress and 
strain in the fiber frame we need to construct a transformation matrix that is used to convert the 
different tensors between the two frames. The components of the transformation matrix are calculated 
using the sine and cosine of   and  . The strain expressed in fiber base is related to the GN strain by 
the following relation [8][9]. 
    
 
2 2
11 11
2 2
22 22
1212
cos sin sin cos
cos sin sin cos
2cos cos 2sin sin sin
d d
d d
dd
     
     
     
     
          
         
  (4) 
Or, we can write  
     d T d    (5) 
Similarly for the stress  
     d T d    (6) 
In the non-orthogonal coordinate system, we can assume that the elastic matrix has an orthotropic 
format. Thus, in this frame, the stress is relied to the strain by:  
  
 
11 21 11
12 21 12 21
12 22 22
12 21 12 21
12
0
1 1
0
1 1
0 0
E E
E E
H
G

   

   

 
  
 
 
  
  
 
 
 
  (7) 
Where  12G  is the in-plane shear rigidity which is a function of shear angle (radians). The shear 
rigidity of the of the glass/polypropylene BPW fabric is given by: 
 
 
4 3 2
12G ( ) 8.48 12.0972 6.1275 0.83 0.051           
     
2.2. Stress computation 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
The aim of this section is to express the 
1n   in the fiber basis. The integration of the equation 
between  
nt  and 1nt   (the time increment  1n nt t t    is assumed to be small enough) using the mid-
point integration scheme gives us: 
 
1n n
  

     (9) 
Where  
 
1/2 1/2
:
n n
H D t
 
     (10) 
Thus, the stresses in the fiber frame are updated according to the Hughes and Winget [4] by the 
following expression  
 
1 1/2 1/2
:
n n n n
H D t 
  
     (11) 
2.3. Implementation of the hypo-elastic model  
The hypo-elastic model described in the previous section has been implemented in a user material 
subroutine VUMAT using FORTRAN code [7]. The implementation is summarized in the following 
figure 2.  According to figure 6, the following matrices and vectors are provided by the code at each 
time step expressed in the default basis of the code. The user subroutine has to return the stress state at 
the end of the step. Then, we compute the inverse matrix of the right stretch tensor. The last tensor is 
used to compute the polar rotation matrix. The current work basis and material basis are computed 
from the tensor F and R. These current vectors are used to calculate the passage matrices T. Then, we 
compute the material strain increment from the code’s strain increment. We update the material stress 
tensor. Finally, we return the stress tensor in the work basis of the code. 
 
Figure 2: Scheme of stress computation in the developed user material subroutine. 
3. Parametric Sensitivity Study 
3.1. Stamping finite element model  
Hemispherical stamping is one of the simplest references of forming woven composites. It has been 
the subject of several publications [11]. The FEM model development of hemispherical and the 
forming tests have been carried out using a commercial finite element code of ABAQUS Explicit 
using standard elements. The model consists of four elements: die, punch, holder and the sheet of 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
composite (figure 3). In this FE model, the three forming tools are considered as rigid solids and are 
modeled as discrete rigid bodies and meshed with 3-D triangular rigid elements (R3D3). The 
reinforcement is considered is considered as a continuum material. Thus, it is modeled as four node 
membrane elements M3D4R.  
 
Figure 3 : Schematic of a deformed plain weave structure with shear deformation 
In the initial configuration, the directions of warp and weft coincide with the two perpendicular edges 
of the membrane element. The material characteristics and process parameters used in this study are 
fixed as [5]: 
- E11=35400MPa  
- E22=35400MPa  
- 12 0    
- Mass density,   = 0.00253 g/mm3 
- Blank holder force F=100N 
- Punch velocity V=90 mm/s 
- Coefficient of friction   = 0.3 
3.2. Effect of variation of the Punch Velocity 
The final state of composite piece is influenced by the process parameters such us velocity of punch, 
binder force, blank holder shape…This section described the effect of the variation of the velocity 
punch on the distribution of the shear angle. In order to control the apparition of many defects like 
plies, wrinkling and buckling, it must select the optimized parameter of punch velocity. The figure 4 
shows a comparison between two punch velocities. 
 
Figure 4: Forming results of shear angle in degrees (SDV71) with punch velocities 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
The lower punch velocity (9mm/s) predicts the shear angle in the final state as 53, which is 1.5 less 
than high speed (90mm/s) shown in Figure 4 
3.3. Effect of variation of Binder force 
The binder or blank holder plays a very important role in the forming processes of composite materials 
like of metallic materials. Mainly, the blank holder is used to ensure a homogeneous flow of the 
material. Thus, it prevents the appearance of manufacturing defects, especially the wrinkling. In order 
to clearly observe the effect of blank holder force on the draped part, two concentrated forces have 
been applied on the binder. The figure 5 shows the results of shear angle with these two different 
forces. The magnitude of the binder force determines correct execution of the forming step. It has a 
little effect on the maximum shear angle for axisymmetric configuration. In fact, the effect of the 
blank holder force is well shown in an asymmetric repartition of force. 
 
 
Figure 5: Forming results of shear angle in degrees (SDV71) with punch velocities 
The figure 6 shows the significant role of the force on the formability results of a composite .In this 
case of configuration, the shear-angle distribution has shown a great variation in the final part. 
 
Figure 6: The evolution and distribution of the shear angle (deg)  in the hemisphere model with 
non-axisymmetric configuration of the holder 
There is also a high value of positive shear angle (73°) along the two free diagonals and (-52°) along 
the two other diagonals  
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
3.4. Effect of variation of the friction coefficient 
Other process parameter which can significantly affect the final geometry is the coefficient of friction 
tools/blank. In this model, there are three interactions: punch/blank, die/blank and holder/blank. These 
interactions are defined by a coefficient of friction introduced in the form of Coulomb friction. In 
reality, the coefficients of friction are different but in this study they are fixed as 0.3. In order to show 
the effect of the friction coefficient over the evolution of the shear angle different values were used. 
Figure shows the sensitivity of this parameter over forming results of two coefficient i.e. 0.3 and 0.5. 
The evolution of positive shear angles is affected by just about 2 but negative shear-angle values differ 
at maximum by almost 3° (Figure 7) 
 
Figure 7: The evolution of the shear angle (degrees) in the two models with  
two different coefficients of friction 
If the friction coefficient is more than 0.5 the results deteriorate largely and are no more 
presentable 
4. Conclusions & outlook 
In this study, numerical simulations of stamping on a dry woven reinforcement have been released 
using a non-orthogonal behavior law which is developed in a non-orthogonal frame (fiber directions) 
in order to keep tracing fiber reorientation during deformation. The hypo elastic model is implemented 
in a commercial FE code (ABAQUS/Explicit) via an user material subroutine VUMAT. In this work, 
we analyses the sensitivity of the shear angle to the process parameters. This parametric study gives 
results that have a good agreement with what is found in the literature. 
Références  
[1] E.R. Fuchs, F.R. Field, R. Roth, R.E. Kirchain, Strategic materials selection in the automobile body: 
economic opportunities for polymer composite design, Compos. Sci.Technol. 68 (9) (2008) 1989–2002. 
[2] S Allaoui, G Hivet, M Haddad, R Agogue, K Khellil, P Beauchene, and Z Aboura. Procédé de mise en 
forme : défauts mésoscopiques et propriétés induites sur le composite. In Journée Scientifique et 
Technique : Défauts dans les composites: Origine, Mesure, Criticité et Impacts sur les Performances. 
[3]: P Badel. Analyse mésoscopique du comportement mécanique des renforts tissés de composites 
utilisant la tomographie aux rayons X. PhD thesis, INSA de Lyon,2008 
[4] Hughes TJR, Winget J. Finite rotation effects in numerical integration of rate constitutive equations 
arising in large deformation analysis. Int J Numer Methods Eng 1980;15:18627. 
[5]CRIESFIELD MA. Nonlinear finite element analysis of solids and structure: advanced topics, vol. 2. 
Chichester: John Wiley Edt.; 1997. 
[6] Muhammad Aurangzeb Khan. Numerical and experimental forming analyses of textile composite 
reinforcements based on a hypo elastic behaviour. 2009. 
[7] Hibbett, Karlsson, and Sorensen. ABAQUS/standard : User's Manual, volume 1. Hibbitt, Karlsson & 
Sorensen, 199 
23
ème
 Congrès Français de Mécanique                              Lille, 28 Août au 1
er
 Septembre 2017 
 
[8] XQ Peng and J Cao. A continuum mechanics-based non-orthogonal constitutive model for woven 
composite fabrics. Composites Part A : Applied Science and Manufacturing, 36(6) :859874, 2005. 
[9] Pu Xue, Xiongqi Peng, and Jian Cao. A non-orthogonal constitutive model for characterizing woven 
composites. Composites part A : Applied Science and manufacturing, 34(2) :183193,2003. 
 


Modélisation et simulation numérique de l'emboutissage d'un renfort tissé sec : Sensibilité de l'angle de cisaillement aux paramètres du procédé

http://documents.irevues.inist.fr/bitstream/2042/63121/1/137400.pdf

Modélisation et simulation numérique de l'emboutissage d'un renfort tissé sec : Sensibilité de l'angle de cisaillement aux paramètres du procédé

Abstract

Similar works

Full text

Available Versions

I-Revues