Analysis and Modeling of Wrinkling in Composite Forming

: Different approaches for the simulation of wrinkling during forming of textile reinforcements are presented. It is shown that 3D ﬁnite element modeling requires the consideration of an additional bending stiffness of the ﬁbers. In shell-type modeling, the bending stiffness is important because it conditions the size of the wrinkles. Different methods to take into account the bending stiffness independently of the tensile stiffness are presented. The onset and development of wrinkles during forming is a global problem that concerns all deformation modes. It is shown in examples that the shear locking angle is not sufﬁcient to conclude about the development of wrinkles. This article highlights the two points common to the different cases of wrinkling of continuous ﬁber textile reinforcements: the quasi-inextensibility of the ﬁbers and the possible slippage between the ﬁbers. It presents and compares different approaches to consider these two aspects. The simulation of the simultaneous forming of multilayered textile reinforcements makes it possible to see the inﬂuence of the orientation of different plies which is an important factor with regard to wrinkling.


Introduction
The forming of textile reinforcements is necessary to make preforms which are then injected in LCM processes (Liquid Composite Molding) [1,2] or to produce composite parts from prepregs [3][4][5]. This forming process conditions the quality of the composite part in particular because it determines the orientation of the fibers. This process must be carried out without defects. Wrinkles are one of the most critical defects in composite manufacturing. It is critical to define the conditions for a forming process that avoids wrinkles [6][7][8][9][10][11][12][13][14][15]. Wrinkles frequently occur when draping the fiber reinforcements. They are possible during the forming of all thin parts, especially metallic ones [16][17][18], but wrinkles are even more frequent for textile reinforcements because their fibrous composition enables slippage between the fibers which greatly reduces bending stiffness. The onset and development of wrinkles were experimentally analyzed in particular by two plane shear tests: the picture frame test and the bias extension test [8,[19][20][21][22][23]. In the picture frame test, the wrinkles appear from a shear angle called shear locking angle which is used as a limit in draping simulations beyond which the wrinkles appear. It is shown in the sequel that the shear angle is not the only quantity that determines wrinkling and that other strains and stiffnesses also play a role. In particular, the bending stiffness determines the size and number of wrinkles.
In order to determine the occurrence of wrinkles and the conditions of a forming process that will prevent these wrinkles, forming simulation codes have been developed over the past three decades [24][25][26][27][28][29] and are regularly improved. They avoid trial-and-error adjustments. The simulation of wrinkles during the forming of a textile reinforcement is currently possible, at least in a large number of cases [7,9,14,30].
It is shown in this paper that 3D simulation models based for example on a hyperelastic model are possible but require a second gradient approach or taking into account the curvature of the fibers to correctly describe the forming of, including the wrinkles. Due to the low thickness of the textile reinforcements of composites, shell-type modeling is the most frequent. It is established that membranes are not suitable for simulating the shape of the wrinkles correctly. In order to consider the bending stiffness which is particular in the case of fibrous reinforcements given the possible slippage between the fibers, several approaches are used, which may be efficient but do not correspond to the physics of the textile material deformation. Within the framework of modeling by stress resultant shells, it is shown that the shear locking angle is insufficient to determine the onset or not of wrinkles. Finally, the analysis of wrinkling during the simultaneous forming of multilayer stacks confirms that a different orientation of the plies of the laminate is a strong source of wrinkles and that this phenomenon can be correctly predicted by simulation. Lastly, the wrinkles created by the bending of a laminate are analyzed.

Kinematic Drape Modeling
The approach is purely geometric. It is based on the following assumptions: There is no slippage between the yarns and the tools after their contact; the yarns constituting the draped fabric are inextensible; the yarns have no bending stiffness; the warp and weft yarns cross each other without slippage and with a free rotation. These methods are also called the fishing net algorithm or pin-jointed net method [31][32][33][34][35]. The method consists of iteratively placing the points of the thread on the tool based on the above assumptions. The determination of a point knowing its neighbors leads to a scalar equation. Consequently, the method is very fast. It is used in some composite design software. Nevertheless, the solution obtained does not consider the material of the textile reinforcement which is draped. On the other hand, the load boundary conditions are not taken into account. The method does not describe wrinkles. The angles between warp and weft that are obtained can give a criterion for the appearance of wrinkles when they exceed a limit angle called shear locking angle. This approach will be discussed in Section 3. Kinematic drape modeling cannot be considered as a method to simulate wrinkles during the forming of a textile reinforcement.

Simulations Based on 3D Elements
To consider the mechanical properties of reinforcements during forming, a 3D finite element modeling can be performed. It is all the more feasible as the textile reinforcement is thick. A specific constitutive model for textile reinforcement during draping is then necessary. Hypoelastic [36][37][38][39][40] and hyperelastic models [41][42][43][44][45] have been developed.
In particular, the hyperelastic model proposed by Charmetant [43] concerns the 3D textile reinforcements. The reinforcement is directed by three-unit vectors M 1 , M 2 , M 3 , respectively, the warp and directions and the direction in the thickness.
Six invariants represent the deformation modes shown in Figure 1: I elongi are the logarithmic strains in the yarn directions M i (Figure 1a,b), I comp is the compaction logarithmic strain (Figure 1c), I cp is the sine of the angle variation between directions M 1 and M 2 (Figure 1d), and I cti are the sine of the angle variation between directions M i and M 3 (Figure 1e,f). The strain energy density is assumed to be defined by the strain energies corresponding to these invariants.
The above invariants are related to the theoretical invariants of the right Cauchy-green tensor C = F T F on which a hyperelastic orthotropic material must depend [43,46]. The strain energy density gives the stress tensor S (Second Piola-Kirchhoff).
. Compos. Sci. 2021, 5, x FOR PEER REVIEW 3 of 17 The above invariants are related to the theoretical invariants of the right Cauchygreen tensor  T C F F on which a hyperelastic orthotropic material must depend [43,46]. The strain energy density gives the stress tensor S (Second Piola-Kirchhoff) Each of the six deformation potentials in Equation (1) is assumed to be polynomial and is identified by tests involving the respective deformation mode [43,47].
3D finite element simulations based on the above hyperelastic behavior can be used to simulate the forming of textile reinforcements, especially when they are thick. Figure 2 shows the hemispherical forming of a self-reinforced polypropylene (SRPP) [48]. The initial SRPP panels were 22 by 22 cm and 2 mm thick. They are obtained from 12 layers of PP twill fabrics stacked alternately with PP films. The hemispherical punch has a diameter of 90 mm. The simulation is performed on a quarter of the panel with symmetry conditions. The other edges are free. Wrinkles develop during forming. They are obtained by simulation using the hyperelastic model described above even if the correspondence of the wrinkling obtained experimentally and by simulation is not perfect. One advantage of the modeling by 3D finite elements is that it considers the deformation in the thickness. A simulation using 3D finite elements of the wrinkles due to the consolidation of prepreg stacks over curved tools is presented in [49]. Each of the six deformation potentials in Equation (1) is assumed to be polynomial and is identified by tests involving the respective deformation mode [43,47].

Displacement
3D finite element simulations based on the above hyperelastic behavior can be used to simulate the forming of textile reinforcements, especially when they are thick. Figure 2 shows the hemispherical forming of a self-reinforced polypropylene (SRPP) [48]. The initial SRPP panels were 22 by 22 cm and 2 mm thick. They are obtained from 12 layers of PP twill fabrics stacked alternately with PP films. The hemispherical punch has a diameter of 90 mm. The simulation is performed on a quarter of the panel with symmetry conditions. The other edges are free. Wrinkles develop during forming. They are obtained by simulation using the hyperelastic model described above even if the correspondence of the wrinkling obtained experimentally and by simulation is not perfect. One advantage of the modeling by 3D finite elements is that it considers the deformation in the thickness. A simulation using 3D finite elements of the wrinkles due to the consolidation of prepreg stacks over curved tools is presented in [49].
The above invariants are related to the theoretical invariants of the right Cauchygreen tensor  T C F F on which a hyperelastic orthotropic material must depend [43,46]. The strain energy density gives the stress tensor S (Second Piola-Kirchhoff) Each of the six deformation potentials in Equation (1) is assumed to be polynomial and is identified by tests involving the respective deformation mode [43,47].
3D finite element simulations based on the above hyperelastic behavior can be used to simulate the forming of textile reinforcements, especially when they are thick. Figure 2 shows the hemispherical forming of a self-reinforced polypropylene (SRPP) [48]. The initial SRPP panels were 22 by 22 cm and 2 mm thick. They are obtained from 12 layers of PP twill fabrics stacked alternately with PP films. The hemispherical punch has a diameter of 90 mm. The simulation is performed on a quarter of the panel with symmetry conditions. The other edges are free. Wrinkles develop during forming. They are obtained by simulation using the hyperelastic model described above even if the correspondence of the wrinkling obtained experimentally and by simulation is not perfect. One advantage of the modeling by 3D finite elements is that it considers the deformation in the thickness. A simulation using 3D finite elements of the wrinkles due to the consolidation of prepreg stacks over curved tools is presented in [49].  Figure 3 shows the hemispherical forming of a thick textile reinforcement. The interlock fabric was 15 × 15 cm and 15 mm thick. The punch has a 100 mm diameter. Figure 3b corresponds to 80% of the simulation. Important wrinkles develop in the preform. These wrinkles are spurious. The forming experiment (Figure 3a) shows no wrinkles. It has been shown that hyperelastic models such as the one presented above, which are a first gradient models, cannot reflect both the very low transverse stiffness of the fibrous material and the bending stiffness of each fiber [50]. To overcome this problem, second gradient approaches have been proposed [51,52]. Another related approach consists in adding a deformation energy corresponding to the bending of the fibers to the deformation energy. This requires the calculation of the fiber curvature which can be performed from the position of the nodes of the neighboring elements ( Figure 3d) [47]. Figure 3c shows the result of the simulation based on this approach which is in good agreement with the experiments and does not show spurious wrinkles.  Figure 3 shows the hemispherical forming of a thick textile reinforcement. The interlock fabric was 15 × 15 cm and 15 mm thick. The punch has a 100 mm diameter. Figure 3b corresponds to 80% of the simulation. Important wrinkles develop in the preform. These wrinkles are spurious. The forming experiment (Figure 3a) shows no wrinkles. It has been shown that hyperelastic models such as the one presented above, which are a first gradient models, cannot reflect both the very low transverse stiffness of the fibrous material and the bending stiffness of each fiber [50]. To overcome this problem, second gradient approaches have been proposed [51,52]. Another related approach consists in adding a deformation energy corresponding to the bending of the fibers to the deformation energy. This requires the calculation of the fiber curvature which can be performed from the position of the nodes of the neighboring elements ( Figure 3d) [47]. Figure 3c shows the result of the simulation based on this approach which is in good agreement with the experiments and does not show spurious wrinkles.

Membrane Approach
Forming simulations based on 3D elements presented above lack numerical efficiency, especially when the thickness of the textile reinforcement is small compared to its dimensions in the plane, which is generally the case. Shell approaches where the dimension in the thickness direction is small are well suited for textile reinforcements of composites. Due to the fibrous composition of the reinforcements, the bending rigidity is low. Membrane approaches that neglect the bending stiffness have been used to simulate the draping of textile fabrics [53][54][55][56][57]. The mechanical behavior in the membrane surface is driven by the warp and weft yarns. For this goal, behavior models sometimes called "nonorthogonal" have been developed [38,53].
Hyperelastic and hypoelastic approaches were also used [36,58,59]. Wrinkles may develop when simulating with membrane elements (Figure 4a). However, these wrinkles are too small and do not represent reality. This is due to the lack of bending stiffness of the membrane approach. The following will show that this bending stiffness, although

Membrane Approach
Forming simulations based on 3D elements presented above lack numerical efficiency, especially when the thickness of the textile reinforcement is small compared to its dimensions in the plane, which is generally the case. Shell approaches where the dimension in the thickness direction is small are well suited for textile reinforcements of composites. Due to the fibrous composition of the reinforcements, the bending rigidity is low. Membrane approaches that neglect the bending stiffness have been used to simulate the draping of textile fabrics [53][54][55][56][57]. The mechanical behavior in the membrane surface is driven by the warp and weft yarns. For this goal, behavior models sometimes called "non-orthogonal" have been developed [38,53].
Hyperelastic and hypoelastic approaches were also used [36,58,59]. Wrinkles may develop when simulating with membrane elements (Figure 4a). However, these wrinkles are too small and do not represent reality. This is due to the lack of bending stiffness of the membrane approach. The following will show that this bending stiffness, although low, is important in the development of the wrinkles and in particular with regard to their size. [9,50]. Bending of the textile reinforcements and the formulation of the corresponding shell approaches is not straightforward. Given the possible slip between the fibers, the deformation is not the same as for standard continuous media. For a given tensile stiffness, the bending stiffness is much lower for textiles than for standard continuous materials. Classical shell theories (Kirchhoff, Mindlin) are not relevant for fibrous media [60]. To solve this difficulty and to have a lower bending stiffness that is decoupled from the tensile stiffness, several approaches have been developed.  Beam finite elements in the direction of the warp and weft yarns can be associated with membrane elements (Figure 5a). The bending stiffness is that of the beams, with the membrane providing the in-plane shear stiffness [61][62][63]. Using finite shell elements and remaining within the framework of classical shell theory, some approaches modify the ratio between bending and tension stiffnesses. The textile reinforcement can be seen as a laminate composed of layers of different Young's modulus so that both tensile and bending stiffnesses are achieved (Figure 5b,c) [64,65]. Figure 5. Introduction of bending stiffness into a textile fabric: (a) combination of beam and membrane elements [62], (b,c) fibrous fabric considered as a laminate [64,65].
The superposition of a membrane finite element and a pure bending shell element is also a way to decouple the two behaviors. In reference [66,67] DKT (Discrete Kirchhoff Triangle) shell element is superimposed on a membrane element Superposition of shell

Decoupling the Membrane and Bending Stiffnesses
The draping shown in Figure 4b was obtained using shell finite elements specific to the fibrous reinforcements. Taking into account the bending stiffness, although low, makes it possible to obtain wrinkles of the appropriate size. To simulate the wrinkles, the bending stiffness must be considered and therefore shell finite elements are required [9,50]. Bending of the textile reinforcements and the formulation of the corresponding shell approaches is not straightforward. Given the possible slip between the fibers, the deformation is not the same as for standard continuous media. For a given tensile stiffness, the bending stiffness is much lower for textiles than for standard continuous materials. Classical shell theories (Kirchhoff, Mindlin) are not relevant for fibrous media [60]. To solve this difficulty and to have a lower bending stiffness that is decoupled from the tensile stiffness, several approaches have been developed.
Beam finite elements in the direction of the warp and weft yarns can be associated with membrane elements (Figure 5a). The bending stiffness is that of the beams, with the membrane providing the in-plane shear stiffness [61][62][63]. Using finite shell elements and remaining within the framework of classical shell theory, some approaches modify the ratio between bending and tension stiffnesses. The textile reinforcement can be seen as a laminate composed of layers of different Young's modulus so that both tensile and bending stiffnesses are achieved (Figure 5b,c) [64,65]. low, is important in the development of the wrinkles and in particular with regard to their size.

Decoupling the Membrane and Bending Stiffnesses
The draping shown in Figure 4b was obtained using shell finite elements specific to the fibrous reinforcements. Taking into account the bending stiffness, although low, makes it possible to obtain wrinkles of the appropriate size. To simulate the wrinkles, the bending stiffness must be considered and therefore shell finite elements are required [9,50]. Bending of the textile reinforcements and the formulation of the corresponding shell approaches is not straightforward. Given the possible slip between the fibers, the deformation is not the same as for standard continuous media. For a given tensile stiffness, the bending stiffness is much lower for textiles than for standard continuous materials. Classical shell theories (Kirchhoff, Mindlin) are not relevant for fibrous media [60]. To solve this difficulty and to have a lower bending stiffness that is decoupled from the tensile stiffness, several approaches have been developed.  Beam finite elements in the direction of the warp and weft yarns can be associated with membrane elements (Figure 5a). The bending stiffness is that of the beams, with the membrane providing the in-plane shear stiffness [61][62][63]. Using finite shell elements and remaining within the framework of classical shell theory, some approaches modify the ratio between bending and tension stiffnesses. The textile reinforcement can be seen as a laminate composed of layers of different Young's modulus so that both tensile and bending stiffnesses are achieved (Figure 5b,c) [64,65]. Figure 5. Introduction of bending stiffness into a textile fabric: (a) combination of beam and membrane elements [62], (b,c) fibrous fabric considered as a laminate [64,65].
The superposition of a membrane finite element and a pure bending shell element is also a way to decouple the two behaviors. In reference [66,67] DKT (Discrete Kirchhoff Triangle) shell element is superimposed on a membrane element Superposition of shell The superposition of a membrane finite element and a pure bending shell element is also a way to decouple the two behaviors. In reference [66,67] DKT (Discrete Kirchhoff Triangle) shell element is superimposed on a membrane element Superposition of shell and membrane elements. This approach enables the modeling of wrinkles during a forming. For example, the wrinkles that develop during the double diaphragm forming process are simulated using this approach in reference [68].
In the approaches described above, the membrane energy and the bending energy are no longer linked and are considered to be decoupled. These methods are nevertheless based on artificial assumptions and are not founded on the physics of textile reinforcement. They provide the deformation within the framework of theories that are not verified for fibrous media.
A shell approach specific to fibrous reinforcements and based on the possible slippage and inextensibility of fibers has been recently proposed in [60,69].

Stress Resultant Shell Elements for Textiles
A stress resultant shell approach for textile reinforcements considers tensions T 11 and T 22 in warp and weft fiber directions, in-plane shear moment C γ and bending moments M 11 and M 22 (Figure 6a). These stress resultants are the conjugates of axial elongation ε 11 , ε 22 , in-plane shear angle γ and curvatures χ 11 , χ 22 . On an RUC (Representative Unit Cell), the internal virtual work is the sum of the virtual works of tension, in-plane shear and bending and takes the following form: ment. They provide the deformation within the framework of theories that are not verified for fibrous media. A shell approach specific to fibrous reinforcements and based on the possible slippage and inextensibility of fibers has been recently proposed in [60,69].

Stress Resultant Shell Elements for Textiles
A stress resultant shell approach for textile reinforcements considers tensions  (Figure 6a). These stress resultants are the conjugates of axial elongation 11 22 ,  , in-plane shear angle  and curvatures 11 22 ,  . On an RUC (Representative Unit Cell), the internal virtual work is the sum of the virtual works of tension, in-plane shear and bending and takes the following form: tension  shear  bending  11  22  11  22  int  int  int  11  1  22  2  11  1 Here 1 L and 2 L are the lengths of the RUC in warp and weft fiber directions. This stress resultant shell approach naturally decouples the virtual work of tension and bending and is therefore well suited to the simulation of the forming of composite textile reinforcements.
(a) (b) The approach requires knowledge of the behaviors in tension, plane shear, and bending, which here are relationships between the stress resultants and the corresponding deformations. Tension tests, picture frame or bias extension test (shear), and bending test enable to determine these behaviors. A triangular finite element with three nodes has been developed in References [70,71]. It uses the position of the nodes of neighboring elements (Figure 6b) to calculate the curvature in the element and take into account the bending deformation energy [72,73]. This shell resultant stress finite element has been shown to be effective for the simulation of the forming of composite textile reinforcements and for the simulation of wrinkles [9,74]. This stress resultant element will be used in the examples presented in the following. Figure 7 shows the experimental and numerical preforms obtained for a stamping of a G1151 interlock fabric (Hexcel) by a tetrahedral punch ( Figure  7a) [75], and the punch has a 28 cm base and a 17.5 cm height. The shear angles obtained Here L 1 and L 2 are the lengths of the RUC in warp and weft fiber directions. This stress resultant shell approach naturally decouples the virtual work of tension and bending and is therefore well suited to the simulation of the forming of composite textile reinforcements.
The approach requires knowledge of the behaviors in tension, plane shear, and bending, which here are relationships between the stress resultants and the corresponding deformations. Tension tests, picture frame or bias extension test (shear), and bending test enable to determine these behaviors. A triangular finite element with three nodes has been developed in References [70,71]. It uses the position of the nodes of neighboring elements (Figure 6b) to calculate the curvature in the element and take into account the bending deformation energy [72,73]. This shell resultant stress finite element has been shown to be effective for the simulation of the forming of composite textile reinforcements and for the simulation of wrinkles [9,74]. This stress resultant element will be used in the examples presented in the following. Figure 7 shows the experimental and numerical preforms obtained for a stamping of a G1151 interlock fabric (Hexcel) by a tetrahedral punch (Figure 7a) [75], and the punch has a 28 cm base and a 17.5 cm height. The shear angles obtained by the simulation are in good agreement with the experimental shear angles (Figure 7b). In particular, the shear angles in the lower corner of the tetrahedral are 60 • (Figure 7c). This value is well obtained by the simulation (Figure 7d). In addition, the wrinkles in the flat part of the preform are simulated in a relevant way. There is no wrinkle in the tetrahedral zone (useful zone) with a blank holder pressure of 1 bar. When this blank holder pressure is decreased to 0.2 bar, the simulation shows that the wrinkles propagate in the tetrahedral zone (Figure 7e). Simulation thus proves to be a means of defining the conditions of the forming process to avoid wrinkles. hedral zone (useful zone) with a blank holder pressure of 1 bar. When this blank holder pressure is decreased to 0.2 bar, the simulation shows that the wrinkles propagate in the tetrahedral zone (Figure 7e). Simulation thus proves to be a means of defining the conditions of the forming process to avoid wrinkles.
This example confirms that wrinkling occurs when the dynamics Equation (4), combined with the boundary conditions, leads to an out-of-plane solution. If the shear deformations are necessary for forming (this is the case for double-curvature geometries) and the shear stiffness becomes large, wrinkling will occur, unless the stresses in the fibers prevent these wrinkling as is the case in Figure 7b in the corner of the tetrahedral.  This example confirms that wrinkling occurs when the dynamics Equation (4), combined with the boundary conditions, leads to an out-of-plane solution. If the shear deformations are necessary for forming (this is the case for double-curvature geometries) and the shear stiffness becomes large, wrinkling will occur, unless the stresses in the fibers prevent these wrinkling as is the case in Figure 7b in the corner of the tetrahedral.

Influence of the Bending Stiffness
The consideration of the bending stiffness is what differentiates the shell approach from the membrane approach. Wrinkles may appear in a membrane model as shown in Figure 4a, but their size is not defined in this approach. The wrinkles are as small as the mesh used allows. If a bending stiffness is added (shell approach), this defines the size of the wrinkles (Figure 4b). Figure 8 shows the simulation of the compression in the fiber direction of a woven sheet with an initial length of 80 to 60 mm. The bending stiffness is set to 1, 5, and 10 Nmm −1 in Figure 8a-c, respectively. The bending stiffness determines the size of the wrinkles and thus their number. One of the roles of the bending stiffness is to determine the size of the wrinkles [9]. The consideration of the bending stiffness is what differentiates the shell approach from the membrane approach. Wrinkles may appear in a membrane model as shown in Figure 4a, but their size is not defined in this approach. The wrinkles are as small as the mesh used allows. If a bending stiffness is added (shell approach), this defines the size of the wrinkles (Figure 4b). Figure 8 shows the simulation of the compression in the fiber direction of a woven sheet with an initial length of 80 to 60 mm. The bending stiffness is set to 1, 5, and 10 Nmm −1 in Figure 8a-c, respectively. The bending stiffness determines the size of the wrinkles and thus their number. One of the roles of the bending stiffness is to determine the size of the wrinkles [9].  Figure 9 shows simulations of the forming of a thermoplastic prepreg [76] in which wrinkles develop. The blank was 40 × 40 cm, and the cylindrical punch has a diameter of 150 mm. The simulations are performed with a bending stiffness which is increased or decreased with respect to the experimental value (Figure 9a). The wrinkles are smaller in size when the stiffness is decreased and larger when the stiffness is increased.  Figure 9 shows simulations of the forming of a thermoplastic prepreg [76] in which wrinkles develop. The blank was 40 × 40 cm, and the cylindrical punch has a diameter of 150 mm. The simulations are performed with a bending stiffness which is increased or decreased with respect to the experimental value (Figure 9a). The wrinkles are smaller in size when the stiffness is decreased and larger when the stiffness is increased. The consideration of the bending stiffness is what differentiates the shell approach from the membrane approach. Wrinkles may appear in a membrane model as shown in Figure 4a, but their size is not defined in this approach. The wrinkles are as small as the mesh used allows. If a bending stiffness is added (shell approach), this defines the size of the wrinkles (Figure 4b). Figure 8 shows the simulation of the compression in the fiber direction of a woven sheet with an initial length of 80 to 60 mm. The bending stiffness is set to 1, 5, and 10 Nmm −1 in Figure 8a-c, respectively. The bending stiffness determines the size of the wrinkles and thus their number. One of the roles of the bending stiffness is to determine the size of the wrinkles [9].  Figure 9 shows simulations of the forming of a thermoplastic prepreg [76] in which wrinkles develop. The blank was 40 × 40 cm, and the cylindrical punch has a diameter of 150 mm. The simulations are performed with a bending stiffness which is increased or decreased with respect to the experimental value (Figure 9a). The wrinkles are smaller in size when the stiffness is decreased and larger when the stiffness is increased. However, it can be seen that the simulations for which the stiffness is divided or multiplied by two (Figure 9b,d) differ little from the reference simulation. On the other hand, the simulation for which the stiffness is divided by ten shows much smaller and more numerous wrinkles (close to the membrane behavior), and the simulation for which the stiffness is multiplied by ten shows much larger and less numerous wrinkles.
To summarize, the bending stiffness is necessary to simulate the wrinkles of the textile reinforcements during forming. The bending stiffness determines the size of the wrinkle. The size of the wrinkles increases with the bending stiffness. However, this dependence is not very sensitive, and a significant change in the bending stiffness is necessary to achieve a significant change in the size and number of wrinkles.

Beyond the Shear Locking Angle
In-plane shear is the main mode of deformation when forming a textile reinforcement. The fibers are very stiff in tension, and elongation is very low. The in-plane shear behavior has been the subject of a large number of studies, in particular with regard to the analysis of the Picture Frame test and the Bias Extension test which have been developed to determine the in-plane shear behavior of textiles [8,[19][20][21][22]77,78]. The in-plane shear stiffness of the textiles is low, which helps draping. Nevertheless, from a certain shear angle, the warp and weft yarns come into contact and the stiffness increases. This increase leads to wrinkles in the picture frame test. The shear angle at which these wrinkles occur is called the "shear locking angle". It is often considered as a characteristic of the textile reinforcement and a value that should not be exceeded during the forming to avoid wrinkles [19,77,[79][80][81]. The increase in in-plane shear stiffness for large shear angles is indeed a favorable factor for wrinkles; however, wrinkling is a more global phenomenon in which all the characteristics of the fibrous reinforcement play a role. This can be seen by considering the virtual work theorem that is written, for any virtual displacement field equal to zero on the boundary with prescribed displacements: in the case of textile reinforcements.
Here, δW ext , δW int , δW acc , are the virtual works of external, internal and acceleration quantities. This equation can lead to in-plane or out-of-plane displacements of the textile fabric, i.e., wrinkles. All the stiffnesses of tension, plane shear, and bending contribute to the solution of Equation (4) and thus to the deformation of the fabric. A large value of the shear angle can lead to wrinkles, but not always. For example, in the lower corner of the tetrahedral in Figure 7, the shear angle is very large (60 • ), but there is no wrinkle in this area. This is due to the effect of the blank holders which create tensions and avoid wrinkles. In other areas of the same textile reinforcement, wrinkles form while the shear angles are much smaller. Boundary conditions are important in the development of wrinkling, especially, the blank holders are intended to prevent wrinkling by imposing tensions in the yarns of the fabric. In Figure 8, the fabric is in compression in the direction of the fibers. There is no in-plane shear. However, the situation leads to significant wrinkles.
Generally speaking, it is the totality of the deformations and stiffnesses in tension, plane shear, and bending that must be considered in the global problem (Equation (4)). A high plane shear stiffness and a state of compression in the direction of the fibers are favorable to the development of wrinkles. A high bending stiffness and a state of tension in the direction of the fibers are favorable for wrinkle-free deformation. For a given geometry and loading, the dynamic Equation (4) will lead or not to wrinkles depending on all the characteristics.
Nevertheless, there are some situations where the comparison of the shear angle (obtained for example by a fishing net algorithm) with the shear locking angle can be used for the prediction of wrinkling. The manual draping of the prepregs can be considered in particular.

Influence of the Orientation of the Plies of a Stack during Forming
Simultaneous forming of a stack of textile reinforcement layers increases the propensity for wrinkle development especially when the plies of the stack have different directions [82][83][84][85][86]. For example, Figure 10 shows  (Figure 10c) leads to significant wrinkles. The latter case is very frequent since it leads to a more isotropic composite. This point is a major difficulty in the forming of multilayered preforms.
Nevertheless, there are some situations where the comparison of the shear angle (obtained for example by a fishing net algorithm) with the shear locking angle can be used for the prediction of wrinkling. The manual draping of the prepregs can be considered in particular.

Influence of the Orientation of the Plies of A Stack during Forming
Simultaneous forming of a stack of textile reinforcement layers increases the propensity for wrinkle development especially when the plies of the stack have different directions [82][83][84][85][86]. For example, Figure 10 shows the hemispherical forming of four-ply stacks of glass plain weave fabric. The fabric stack is 30 × 30 cm, and the diameter of the punch is 150 mm. Stack forming where the wrinkles have the same orientation, [0°/90°]4 ( Figure  10a) or [−45°/45°]4 (Figure 10b), is performed without wrinkles. On the other hand, stack forming which is composed of plies of different orientations [0°/90°,−45°/45°]2 (Figure 10c) leads to significant wrinkles. The latter case is very frequent since it leads to a more isotropic composite. This point is a major difficulty in the forming of multilayered preforms.
A simulation of these multilayer forming processes based on the stress resultant shell presented in Section 2.4.2 leads to wrinkles or not in accordance with the experiments (Figure 11). Simulation can be used to determine if the forming will be wrinkle free or not for given orientations of the stack.  A simulation of these multilayer forming processes based on the stress resultant shell presented in Section 2.4.2 leads to wrinkles or not in accordance with the experiments (Figure 11). Simulation can be used to determine if the forming will be wrinkle free or not for given orientations of the stack. To understand the formation of wrinkles when the layers are oriented differently, the value of the tensions in a 45° oriented ply are analyzed. Figure 12a shows the tensions in the warp direction which are all positive (traction). On the other hand, a part of the tensions in the weft direction (Figure 12b) is negative (compression). These compressions in the weft yarn direction lead to wrinkles [86]. They are due to the relative slippage with Figure 11. Numerical deformed geometry of four plies formed simultaneously for a stack with all the plies with the same orientation and for a ply with different orientation.
To understand the formation of wrinkles when the layers are oriented differently, the value of the tensions in a 45 • oriented ply are analyzed. Figure 12a shows the tensions in the warp direction which are all positive (traction). On the other hand, a part of the tensions in the weft direction ( Figure 12b) is negative (compression). These compressions in the weft yarn direction lead to wrinkles [86]. They are due to the relative slippage with friction between the layer and its neighbors that are necessary for forming due to the different orientation of the plies. These compressions in the direction of the weft yarn are related to the friction between the layers. Figure 13 shows the influence of the coefficient of friction on the amplitude of these compressions. When the coefficient of friction is zero, there is no compression in the weft yarn (Figure 13a). Compressions appear with a coefficient of friction of 0.2 ( Figure 13b) and increase for a coefficient of friction of 0.5 (Figure 13c). Decreasing the coefficient of friction between layers is not easy. Nevertheless, an active system has been proposed to control the separation of the layers during forming and thus avoid friction [87]. Figure 11. Numerical deformed geometry of four plies formed simultaneously for a stack with all the plies with the same orientation and for a ply with different orientation.
To understand the formation of wrinkles when the layers are oriented differently, the value of the tensions in a 45° oriented ply are analyzed. Figure 12a shows the tensions in the warp direction which are all positive (traction). On the other hand, a part of the tensions in the weft direction ( Figure 12b) is negative (compression). These compressions in the weft yarn direction lead to wrinkles [86]. They are due to the relative slippage with friction between the layer and its neighbors that are necessary for forming due to the different orientation of the plies. These compressions in the direction of the weft yarn are related to the friction between the layers. Figure 13 shows the influence of the coefficient of friction on the amplitude of these compressions. When the coefficient of friction is zero, there is no compression in the weft yarn (Figure 13a). Compressions appear with a coefficient of friction of 0.2 ( Figure 13b) and increase for a coefficient of friction of 0.5 ( Figure  13c). Decreasing the coefficient of friction between layers is not easy. Nevertheless, an active system has been proposed to control the separation of the layers during forming and thus avoid friction [87].

Analysis of Wrinkling during Bending of Multilayer Composite Reinforcements
The manufacturing processes for preforms and prepregs often involve folding or bending of laminates. Depending on the boundary conditions, these bending operations are conducive to the development of wrinkles in the laminate [49,88,89]. It is important to be able to predict the onset and development of these wrinkles through numerical simulations. The stress resultant shell finite elements presented in Section 2.4.2 prove their efficiency in the analysis of wrinkles induced by bending.

Analysis of Wrinkling during Bending of Multilayer Composite Reinforcements
The manufacturing processes for preforms and prepregs often involve folding or bending of laminates. Depending on the boundary conditions, these bending operations are conducive to the development of wrinkles in the laminate [49,88,89]. It is important to be able to predict the onset and development of these wrinkles through numerical simulations. The stress resultant shell finite elements presented in Section 2.4.2 prove their efficiency in the analysis of wrinkles induced by bending.

Influence of the Stack Sequence
A multilayer composite reinforcement was bended by a 45 • rotation of its two edges which are clamped in the support (Figure 14). The lay-up sequence of different layers of multilayer composite reinforcements has a significant influence on the development of wrinkles. Figure 10 shows the experiments and simulations of the deformed shapes of multilayer stack under different orientation of textile reinforcements during the symmetrical bending [89]. As the layers are oriented in [0 • /90 • ] 10 , the inextensibility of fibrous reinforcement results in a large wrinkle during the bending process (Figure 14a). When the layers of the stack were all oriented at [−45 • /+45 • ] 10 , a wrinkle-free deformed shape of the laminate was obtained due to the in-plane shear of the reinforcement while bending (Figure 14b). In the case where the [0 • /90 • ] and [−45 • /+45 • ] plies are mutually crossed, the morphology of the wrinkle is changed but the presence of [0 • /90 • ] layer makes the wrinkle still exist, but one out of every two plies (Figure 14c). Simulations of the bending of the stacks well describe these wrinkle phenomena in different cases and can be used to optimize processes and stacks to avoid wrinkles during bending operations.

Analysis of Wrinkling during Bending of Multilayer Composite Reinforcements
The manufacturing processes for preforms and prepregs often involve folding or bending of laminates. Depending on the boundary conditions, these bending operations are conducive to the development of wrinkles in the laminate [49,88,89]. It is important to be able to predict the onset and development of these wrinkles through numerical simulations. The stress resultant shell finite elements presented in Section 2.4.2 prove their efficiency in the analysis of wrinkles induced by bending.

Influence of the Stack Sequence
A multilayer composite reinforcement was bended by a 45° rotation of its two edges which are clamped in the support (Figure 14). The lay-up sequence of different layers of multilayer composite reinforcements has a significant influence on the development of wrinkles. Figure 10 shows the experiments and simulations of the deformed shapes of multilayer stack under different orientation of textile reinforcements during the symmetrical bending [89]. As the layers are oriented in [0°/90°]10, the inextensibility of fibrous reinforcement results in a large wrinkle during the bending process (Figure 14a). When the layers of the stack were all oriented at [−45°/+45°]10, a wrinkle-free deformed shape of the laminate was obtained due to the in-plane shear of the reinforcement while bending (Figure 14b). In the case where the [0°/90°] and [−45°/+45°] plies are mutually crossed, the morphology of the wrinkle is changed but the presence of [0°/90°] layer makes the wrinkle still exist, but one out of every two plies (Figure 14c). Simulations of the bending of the stacks well describe these wrinkle phenomena in different cases and can be used to optimize processes and stacks to avoid wrinkles during bending operations.

Influence of the Out-of-Plane Pressure
The generation of wrinkles during the bending of multilayer composite reinforcements is strongly associated with the out-of-plane pressure to which they are subjected. Figure 15 presents the deformation of multilayer stack during the L-flange forming under different loads [89][90][91][92]. When a small load (1 N) was applied at the left end of the reinforcements, no wrinkles were created, and a natural bevel was formed (Figure 15a).

Influence of the Out-of-Plane Pressure
The generation of wrinkles during the bending of multilayer composite reinforcements is strongly associated with the out-of-plane pressure to which they are subjected. Figure 15 presents the deformation of multilayer stack during the L-flange forming under different loads [89][90][91][92]. When a small load (1 N) was applied at the left end of the reinforcements, no wrinkles were created, and a natural bevel was formed (Figure 15a). As the compression load at the left end continued to increase from 0 N to 20 N, the slippage between layers was constrained, leading to the generation of wrinkles during the bending process (Figure 15b,c). These phenomena are well depicted in the simulation by the stress resultant shell elements presented in Section 2.4.2. Other cases of analysis of wrinkles during laminate bending are given in [89]. As the compression load at the left end continued to increase from 0 N to 20 N, the slippage between layers was constrained, leading to the generation of wrinkles during the bending process (Figure 15b,c). These phenomena are well depicted in the simulation by the stress resultant shell elements presented in Section 2.4.2. Other cases of analysis of wrinkles during laminate bending are given in [89].

Conclusions
Wrinkles are one of the main defects that can occur during the forming of composite preforms. The low bending stiffness of textiles is favorable to wrinkling. It is important that the simulations accurately describe the onsets and development of wrinkles. Shell-type modeling for textile reinforcements must consider the bending stiffness of fibrous media which is very particular given the possible slippage between fibers. The classical theories of Kirchhoff and Mindlin are not suitable for fiber reinforcements and cannot be used. The bending stiffness plays an important role since it determines the size and number of wrinkles (among other things). Several approaches have been proposed to decouple bending and membrane stiffnesses. They remain quite far from the deformation physics of textile reinforcements, the two main aspects of which are the quasi-inextensibility of the fibers and the possible slippage between the fibers. The stress resultant shell approach is more natural. It allows to separate the tension and bending behaviors and is efficient to simulate the onset and development of wrinkling of textile reinforcements. However, it does not give the distribution of deformations and stresses in the thickness of the reinforcement. A shell modeling specific to fibrous media should provide an answer to this aspect.
All strain energies (tension, in-plane shear, bending) play a role in the global dynamic equation and in the formation of the wrinkles. In particular the shear angle and its comparison to the shear locking angle is not sufficient in forming processes.
In the case of multilayered reinforcements, the orientation of the layers plays an important role in the development of wrinkling. When these orientations are different, significant slippage between the layers is created by the forming process, leading to compressive forces in the yarns and wrinkling. High friction between the layers increases this phenomenon. In all of these wrinkling situations, simulation provides major information. In particular, it makes it possible to determine forces and stresses in the fibers of the different layers of the reinforcement that are very difficult (or impossible) to measure.