# Numerical Tools for Composite Woven Fabric Preforming

Authors: Abel Cherouat; Houman Bourouchaki
Journal: Advances in Materials Science and Engineering
 (2013)

---

## Abstract

An important step in the manufacturing processes of thin composite components is the layingup of the reinforcement onto the mould surface. The prediction of the angular distortion of the reinforcement during draping and the changes in fibre orientation are essential for the understanding of the manufacture process and the evaluation of the mechanical properties of the composite structures. This paper presents an optimization-based method for the simulation of the forming processes of woven fabric reinforced composites. Two different approaches are proposed for the simulation of the draping of woven fabric onto complex geometries: geometrical and mechanical approaches. The geometrical approach is based on a fishnet model. It is well adapted to predimensioning fabrics and to give a suitable quantification of the resulting flat patterns. The mechanical approach is based on a mesostructural model. It allows us to take into account the mechanical properties of fibres and resin and the various dominating mode of deformation of woven fabrics during the forming process. Some numerical simulations of the forming process are proposed and compared with the experimental results in order to demonstrate the efficiency of our approaches.


---

## Body

## 1. Introduction

Composite reinforced by woven fabric is known to have high specific stiffness and, in combination with automatic manufacturing processes, makes it possible to fabricate complex components (aircraft, boat, automotive, and military) with high level of weight and cost efficiency due to their excellent mechanical properties. One of the most attractive properties of the fabric reinforcements is that they are easily handled and automatically processed, which can lower the production cost. As known, the substitution of metal alloys by composite materials, in general, reduces structural mass by 20–30%. The mass increase is also due to the numerous variety of semiproducts (roving, fabrics, knitted fabrics, and braids preimpregnated or not) permitting the development of new structures. Fabrication processes, also, have undergone substantial evolution in recent years. Although the traditional layup process will remain the process of choice for some applications, new developments in resin transfer molding (RTM) or sheet molding compound (SMC), low temperature curing prepregs, and low pressure molding compounds have matured significantly are reached, and are now being exploited in high technology areas such as aerospace industry. For example, by using such composites, the automotive industry can realize improved fuel economy through vehicle-weight reduction by replacing the currently used steel and aluminum parts with thermoplastic composites with the added benefit of a corrosion-resistant material. The choice of the manufacturing process depends on the type of matrix and fabric, the temperature required to form the part, and the cost effectiveness of the process. In particular, thermoforming is a promising manufacturing process for producing high-volume low-cost composite parts using commingled fibre glass/polypropylene woven fabrics [1].The simulation of the manufacturing of a textile-reinforced composite part with a liquid composite molding-like process, which involves draping (or deep drawing) and impregnation of the preform, includes several stages (see Figure1). First, a mould is designed with CAD software, and the CAD model is meshed. Then, a draping (or deep drawing) simulation tool is used to compute the deformations of the textile layers inside the mould. As a result, for every element of the mesh, textile parameters like the shear angle and the thickness of the layer are available. With these parameters given, the local (mesoscale) properties for every element are determined (preprocessing) and given as input for the macroscopic structure after resin polymerization simulation. The result of the macrosimulation is then postprocessed to optimize the mechanical properties of the composite structure. The numerical simulation of composite forming is an efficient means of evaluating the factors related to the manufacturing processes and an efficient help to design preforming sequence for the manufacturing of fabric-reinforced composites. It is possible to detect the main problem occurring during the shaping deformation and to obtain good quantitative information on the forming process [2–5].Figure 1
Illustration of the manufacturing of a textile reinforced composite part.Different levels of modeling intervene in the simulation of woven reinforcement forming: (1) architecture design level, (2) preliminary design level, (3) mechanical level by computational software, and (4) optimization level [4–6]. Most of these levels are integrated and take into account specific constraints of the manufacturing processes. The particular form of composite fabrication (prepreg) begins with the preimpregnation of reinforcement materials with a resin. The combining of these two materials occurs prior to the molding process and therefore enables a very accurate reinforcement to resin ratio to be achieved. Prepreg materials are used extensively in the aerospace or automotive industry due to their ability to maximize strength to weight ratios. Prepregs are pliable and therefore are able to be cut into various shapes or patterns prior to processing into the molded products. But for the manufacturing of nondevelopable composite part (part that cannot be made flat and unstretched), a new problem intervenes in the design chain resulting from the number of parameters influencing the global behavior of composite forming process. The ability to define, in advance, the ply shapes and material orientation allowed the engineers to optimize the composite structural properties of the composite products for maximum strength, maximum material utilization, and maximum layup efficiency [7–15].The composite manufacturing process involves large displacements and rotations and large shear of weft and warp fibres, which can have a significant effect on the processing and structural properties of the finished product. The effective modeling of the forming properties of woven fabric polymer composite materials has been the focus of research for more than a decade. The formulation of new and more efficient numerical models for the simulation of the shaping composite processes must allow for reduction in the delay in manufacturing of complex parts and an optimization of costs in an integrated design approach [15–20]. Several modeling approaches have been developed to account for the evolution of the orthotropic directions during high shearing, and these approaches include the geometrical and the finite element approaches.In this context, simulation methods are needed to anticipate the performance of the final part but also to predict the reinforcement preforming and the resin injection. Early methods, based on a geometric approach of the fishnet type algorithm, emerged a few decades ago. The geometrical approach based on the so-called fishnet algorithms is used to determine the deformed shape of draped fabrics. The main assumptions are that the fibers are inextensible, the intersection points between warp and weft yarns are fixed during preforming, and the rotations between warp and weft yarns are free. This method, where the fabric is placed progressively from an initial line, provides a close enough resemblance to handmade draping [18–21]. These methods were originally developed for prepreg hand draping. They are very fast and fairly efficient in many prepreg draping cases. Nevertheless, this method has major drawbacks. They account neither for the mechanical behavior of the fabric nor for the static boundary conditions. This last point is very important in the case of forming with punch and die (such as in the preforming of the RTM process). The loads on the tools, especially on the blank holder, influence the quality of the shaping operation and therefore need to be considered in simulations [21, 22]. More recently finite element methods have been used for the draping simulation of composite reinforcements.The alternative to the geometrical approach consists of a mechanical analysis of the fabric deformation under the boundary conditions prescribed by the forming process. This requires a specific model of the woven reinforcement and its mechanical behavior. The mechanical behavior of woven fabrics is complex due to the intricate interactions of the fibres. It is a multiscale problem. The macroscopic behavior is very much dependent on the interactions of fibres at the mesoscale (scale of the woven unit cell) and at the microscale (level of the fibres constituting yarns). Despite of the great amount of work in the field, there is no widely accepted model that accurately describes all the main aspects of fabric mechanical behavior. The main model families come from the multiscale nature of the textile. A first family of models is obtained by homogenizing the mechanical behavior of the underlying mesostructure and considering the fabric as an anisotropic continuum [9–15]. If these models can easily be integrated in standard finite element using conventional shell or membrane elements, then the identification of homogenized material parameters is difficult, especially because these parameters change when the fabric is strained and when, consequently, the directions and the geometry (crimp, transverse sections, etc.) of the fibres change. Some of these approaches will be described, especially a nonorthogonal constitutive model [23–31] and an anisotropic hypoelastic continuous behavior for fibrous material based on an objective derivative using the rotation of the fibre [24]. Conversely, some authors present finite element models to model the warp and the weft fibres behavior [28]. Each yarn or each fibre is modeled and is assumed to be a straight or a curved beam or truss. Sometimes they are modeled as 3D domains [28]. Springs are often used to model warp and weft yarn interactions. In the objective of fabric forming simulations, some authors extend the finite element modeling to the whole textile structure that is represented by a network of interwoven trusses or beams with different tensional and rotational springs. Accounting for the simplicity of each component, the whole textile structure deformation can be computed.Nevertheless, the computational effort needed is relatively significant. At present, this method is restricted to simple geometry of the local yarn and relatively simple mechanical behavior. When a fine model of the fibrous yarns is used, the analysis can only consider a small part of the textile reinforcement such as a few woven or knitted cells. The semidiscrete approach is a compromise between the above continuous and discrete approaches [14]. A finite element method is associated to a mesoscopic analysis of the woven unit cell. Specific finite elements are defined and are made of a discrete number of woven unit cells. The mechanical behaviour of these woven cells is obtained by experimental analyses or from 3D FE computations of the woven cell. The nodal interior loads are deduced from this local behavior and the corresponding strain energy in the element deformation.For experimental technique, Kenkare and May-Plumlee  [32] distinguished draping from paperiness (paper bending behaviour) by accounting for 3D deformation behavior using a drapemeter that is capable of distorting fabric samples in all three dimensions over a circular disc. Providing an accurate vertical projection of the draped sample by optical means allowed the determination of the drape coefficient (DC), which they defined as “the fraction of the area of the annular ring covered by the projection of the draped sample.” Hence, a high drape coefficient reflects poor or little deformation.Over several years, it was the main concern of drape researchers to accurately record the draped contour. From complex projection techniques [33] to the use of photovoltaic cells [34] and “cut-and-weigh” techniques, image analysis techniques [35] have definitely simplified investigations, providing higher accuracies. Most measuring methods produced similar results. However, image analysis techniques greatly reduced result dispersion as well as measuring time.In this study, geometrical and finite element approaches are used to simulate the deformation of preforming of woven fabric. These approaches, while giving good results and being efficient in terms of computing time, are generally somewhat complex and sometimes very challenging to implement into commercially available FEA packages.(i)
The geometrical approach is well adapted to preliminary design level. It is based on geometrical aspects of the warping. It is based on a modified MOSAIC algorithm, which is suitable to generate a regular quad mesh representing the layup of the curved surfaces (giving the exact fibre orientations). The method is implemented in the GeomDrap software which is now integrated in the ESI-Pam software. This software provides a fibre quality chart (showing the fibre distortions, the rate of falling, and the rate of draped surface) to predict local folding due to that the overlapping of fibers in the shear exceeds limit value (up to 60° in some cases). It can be used to optimize the draping process (with respect to the above quality measure) by improving the layup directions or the marker data location. The layup of complex curved surfaces can be made in a few seconds [18, 22].(ii)
The mechanical approach is based on a mesostructural description for finite deformations and geometrical nonlinearity. The unpolymerized resin has a viscous behavior, and the reinforced fibers are treated as either unidirectional or nonlinear elastic behavior. The unit cell of the mesoscopic model used here for a plain weave fabric consists of bicomponent finite elements. The tensile load is carried by the 1D elements that will capture the changes in the orthotropic directions during the shearing. The 2D element accounts only for the shearing resistance of the fabric and hence has no tensile stiffness. The bicomponent finite elements for modeling composite fabric behavior are based on 3D membrane finite elements representative of resin behavior and truss finite elements representative of warp and weft fibers behavior. The efficiency of the proposed model resides in the simplicity of its finite element discretization and the performance of its mechanical background [30, 31].


## 2. Geometrical Approach

The draping of composite fabric using a mechanical approach requires the resolution of equilibrium PDE’s problems by the finite element method (see Section3). In general, in the case of complex surfaces, the boundary conditions are not well defined, and the contact between the surface and the fabric is difficult to manage. Furthermore, the resolution of such a problem can be too long in CPU time and is detrimental to the optimization stage of draping regarding the initial fibre directions. All of these facts lead us to consider rather a geometrical approach which is very fast and more robust allowing simultaneously to define the stratification sequences and the flat pattern for different plies and to predict difficult impregnated areas which involve manual operations like dart insertion or, on the contrary, the shortage of fabric. Based on technical criteria (mold surface covering, fabric drape covering, and fibre angular distortion), this approach can constitute the predimensioning or the preoptimization stage for the manufacturing of complex composite parts. The geometrical approach is based, in general, on the fishnet method for which a fabric mesh element is subjected only to shear deformations. The difficulty of such a method is the mapping of the fabric mesh element onto any surface. Within this context, several algorithms approximating the geometry of a fabric mesh element plotted onto the surface are proposed. In particular, the edges of the fabric mesh element are approximated by line segments representing a pure estimate in the highly curved area (see [19–22] for a synthesis).In this paper, we propose a new geometrical algorithm which takes into account the true geometry of the fabric mesh element plotted onto the surface. Such a fabric mesh element is then defined by a curved quadrilateral whose edges are geodesic lines with the same length plotted onto the surface to drape. Given three vertices of the fabric mesh element on the surface, we propose an optimization algorithm to define the fourth vertex of the fabric mesh element. This algorithm allows us to drape the complex surface using an advancing front approach from the data of an initial start point between the fabric and the surface and the initial fibre directions at this point. In this section, geometrical approach based on the mould mesh discretization is presented. For the second formulation, we propose an algorithm of composite fabric draping without any approximation on the geometry of surface to be draped.


### 2.1. Geometrical Formulation

First, we present the mathematical formulation of the geometrical draping, and then we propose an algorithm scheme to solve the draping problem. LetΣ denote the surface of the part to drape, and we assume that a geometrical mesh TΣ of surface is known. Let Φ be the woven composite fabric modeled by two families (warp and weft) of mutually orthogonal and inextensible fibre described by the local coordinates x=(ξ,η). These families constitute regular quadrilateral fabric mesh TF of the fabric Φ (Figure 2 gives example of draping steps of complex part). The problem of geometrical draping of Φ onto the surface Σ consists of calculating each node displacement of fabric mesh TF with a point of the surface mesh TΣ such that the lengths of the edge of the corresponding mesh TFΣ on the surface are preserved (not extensible). This problem presents infinity of solutions depending on(i)
starting point associated with a node of fabricTFΣ,(ii)
initial warp and weft orientationα.Geometrical draping steps of woven composite fabric.
(a)
Starting point and draping direction
(b)
Geodesic lines and fibres directions
(c)
Computing the four nodes of fabric cell
(d)
Flat pattern of the initial fabric
(e)
Shear angles between fibres
(f)
Stretching of fibresThus, to ensure a unique solution, we suppose that the points of impact on the part surface as well as the fabric orientation are given. The draping scheme is given by the following step [36]:(1)
associate a starting point (corresponding to the point of impact of the machine to drape) on the surface on geometrical part meshx0Σ=(ξ0,η0) (Figure 2(a)),(2)
compute numerically step by step the warp nodes ofTFΣ, classified as α-nodes, from the starting point, associated with nodes (ξ,η0) of TF (Figure 2(b)),(3)
compute numerically step by step the weft nodes ofTFΣ, classified also as α-nodes, from the starting point, associated with nodes (ξ0,η) of TF (Figure 2(c)),(4)
compute numerically cell by cell all the other nodes ofTFΣ, classified as β-nodes, from x0 and the nodes associated with nodes (ξ,η0) and (ξ0,η) of TF (Figure 2(c)).The nodes ofTFΣ associated with nodes (ξ,η0) and (ξ0,η) of TF and the α-nodes are located on the surface along the geodesic lines emanating from the point of impact. Regarding the β-nodes, various algorithms are proposed [20]. Most of them use an analytical expression of the surface and formulate the draping problem in terms of nonlinear partial differential equations. Other algorithms are also proposed to simplify these equations by using a finite element discretization of the surface by flat triangular face (i.e., a mesh of the surface). Based on this latter approach, we propose a new algorithm. The β-nodes are computed by solving an optimization problem corresponding to determine a vertex of an equilateral quad plotted on the surface from the data of the three other vertices. This optimization problem formulates the direction of the geodesic lines emanating from the searched vertex.Consequently, two problems arise:(1)
problem 1: determine the geodetic exit of a given point of surface according to a given orientation,(2)
problem 2: determine the geodetic exits of these points intersecting themselves mutually according to given two points of surface and lengths (this geodetic is given according to their orientations).The developed algorithm is implemented in GeomDrape tool [18]. This software provides a fibre quality chart (showing the fibre distortions, the rate of falling, and the rate of draped surface) to predict difficult impregnated regions.  It can be used to optimise the draping process (with respect to the above quality measure) by improving the layup directions or the marker data location.  The layup of complex curved surfaces can be made in a few seconds. The use of draping tool should not only allow for a more detailed assessment of the draped fabric (Figure 2(c)), but also for an optimization in terms of total fabric shear deformation (Figure 2(e)), flat pattern of the initial fabric (Figure 2(d)) and fibre stretching (Figure 2(f)), and the optimization of the best layup start point (Figure 2(a)).


### 2.2. Application of Geometrical Approach

Two draping simulation examples are given. These simulations are performed using the geometrical analysis computer code GeomDrap [18]. For each example, we assume that a mesh of the piece to drape is given. The first example shows the influence of the fibre orientations draping in the draping process. The second example shows the efficiency of the proposed method to simulate geometrically the draping composite part.The first example concerns the geometrical draping of a complex piece composed of two half hemispheres with a radius of 38.8 mm that are connected by a half cylinder with a length of 170 mm [23]. Figure 3 shows the piece as well as a mesh of the mechanical part. The centroid of the part is chosen as the impact point for which two different fibre orientations are specified. Figure 4 shows the resulting 3D surface layup for the 0°/90° fibre orientation (Figure 4(a)) and the corresponding 2D flat patterns (Figure 4(b)). Likewise, Figure 5 shows the draping results for the ±45° fibre orientation (Figure 5(a) 3D surface layup and Figure 5(b) the corresponding 2D flat pattern). The shears angles between fibres are given in Figure 6. One can notice that, in the considered cases, the surface of the piece is draped globally. However, to drape completely the proposed surface without excessive fibre distortion <50°, it is necessary to make cuts chisel along the direction 45° for 0°/90° fibre orientation (Figure 6(a)) and along the direction 0°  ± 45° fibre orientation (Figure 6(b)). In case we can allow excessive distortions of the fibers without cutting chisel, the results given in Figure 7 show the distribution of fibers in both directions 0°/90° (Figure 7(a)) and ±45° (Figure 7(a)). We can see that the shear limit reached 52° in the case of 0°/90° fibre orientation and 89° in the case of ±45° fibre orientation.Figure 3
CAD of complex part.Draping of (0°/90°) fabric: (a) 3D surface layup and (b) corresponding 2D flat pattern.
(a)
(b)Draping of (−45°/45°) fabric: (a) 3D surface layup and (b) corresponding 2D flat pattern.
(a)
(b)Shear angles: (a) 0°/90° fabric orientation and (b) −45°/+45° fabric orientation.
(a)
(b)Optimal 3D surface layup: (a) 0°/90° fabric orientation and (b) −45°/+45° fabric orientation.
(a)
(b)To compare the experimentally obtained fibre angle result [23] with the geometrical predicted ones, two cross sections where most shearing occurs are examined. The amount of shear along the edge of these cross sections is used to compare with the draping approach. Figure 3 gives an overview of the two cross sections along the symmetrical L1 and along the median line L2 for different ply orientations. The line L2 corresponds with the examined cross section for a 0°/90° fabric draping orientation and the line L1 for a 45 ± 45° fabric draping orientation. From Figure 8(a), it can be concluded that for the 0°/90° orientation, the agreement between the experimental shear angles and the geometrical draping results is good. On the other side, the drape orientation ±45°, the results from the geometrical model do not agree at all with the experimental results [23] (see Figure 8(b)). The oversimplification of the fabric deformation in this geometrical model gives shear angles up to 89° in case of the 45° fabric orientation, which is impossible for fabric composite. The geometrical model is used with a cutoff shear angle based either on an experimentally determined locking angle or the maximum orientation that the designer is prepared to tolerate, when defining a cutoff angle of 38°, which equals the experimentally determined locking angle. From the above discussed results, it can be concluded that the geometrical draping method is not appropriate for predicting the fiber directions of a fabric sheet after stamp forming.Shear angles of optimal geometrical draping: (a) along the line L2 of 0°/90° and (b) along the line L1 of −45°/45°.
(a)
(b)The second example is the draping of a complex shape (car hood). The centroid of the part is chosen as the starting point from which the 0°/90° and ±45° fibre orientations are specified. Figure9 shows the resulting 3D draping for the two orientations (Figure 9(a) for 0°/90° fabric orientation and Figure 9(b) for 45°/45° fabric orientation). We can note that all the part surface is completely draped. Figure 10 presents shaded contours interpolated from the map of the fiber distortions of 0°/90° (Figure 10(a)) and ±45° fiber orientations (Figure 10(b)). The fiber distortions for both 0°/90° and ±45° draping are very small but the maximum shear angle localization is different.3D geometrical draping of hood car part: (a) 0°/90° fabric orientation and (b) −45°/45° fabric orientation.
(a)
(0°/90°)
(b)
(−45°/45°)Isovalues of shear angles between fibres: (a) 0°/90° fabric orientation and (b) −45°/45° fabric orientation.
(a)
(0°/90°)
(b)
(−45°/45°)


## 2.1. Geometrical Formulation

First, we present the mathematical formulation of the geometrical draping, and then we propose an algorithm scheme to solve the draping problem. LetΣ denote the surface of the part to drape, and we assume that a geometrical mesh TΣ of surface is known. Let Φ be the woven composite fabric modeled by two families (warp and weft) of mutually orthogonal and inextensible fibre described by the local coordinates x=(ξ,η). These families constitute regular quadrilateral fabric mesh TF of the fabric Φ (Figure 2 gives example of draping steps of complex part). The problem of geometrical draping of Φ onto the surface Σ consists of calculating each node displacement of fabric mesh TF with a point of the surface mesh TΣ such that the lengths of the edge of the corresponding mesh TFΣ on the surface are preserved (not extensible). This problem presents infinity of solutions depending on(i)
starting point associated with a node of fabricTFΣ,(ii)
initial warp and weft orientationα.Geometrical draping steps of woven composite fabric.
(a)
Starting point and draping direction
(b)
Geodesic lines and fibres directions
(c)
Computing the four nodes of fabric cell
(d)
Flat pattern of the initial fabric
(e)
Shear angles between fibres
(f)
Stretching of fibresThus, to ensure a unique solution, we suppose that the points of impact on the part surface as well as the fabric orientation are given. The draping scheme is given by the following step [36]:(1)
associate a starting point (corresponding to the point of impact of the machine to drape) on the surface on geometrical part meshx0Σ=(ξ0,η0) (Figure 2(a)),(2)
compute numerically step by step the warp nodes ofTFΣ, classified as α-nodes, from the starting point, associated with nodes (ξ,η0) of TF (Figure 2(b)),(3)
compute numerically step by step the weft nodes ofTFΣ, classified also as α-nodes, from the starting point, associated with nodes (ξ0,η) of TF (Figure 2(c)),(4)
compute numerically cell by cell all the other nodes ofTFΣ, classified as β-nodes, from x0 and the nodes associated with nodes (ξ,η0) and (ξ0,η) of TF (Figure 2(c)).The nodes ofTFΣ associated with nodes (ξ,η0) and (ξ0,η) of TF and the α-nodes are located on the surface along the geodesic lines emanating from the point of impact. Regarding the β-nodes, various algorithms are proposed [20]. Most of them use an analytical expression of the surface and formulate the draping problem in terms of nonlinear partial differential equations. Other algorithms are also proposed to simplify these equations by using a finite element discretization of the surface by flat triangular face (i.e., a mesh of the surface). Based on this latter approach, we propose a new algorithm. The β-nodes are computed by solving an optimization problem corresponding to determine a vertex of an equilateral quad plotted on the surface from the data of the three other vertices. This optimization problem formulates the direction of the geodesic lines emanating from the searched vertex.Consequently, two problems arise:(1)
problem 1: determine the geodetic exit of a given point of surface according to a given orientation,(2)
problem 2: determine the geodetic exits of these points intersecting themselves mutually according to given two points of surface and lengths (this geodetic is given according to their orientations).The developed algorithm is implemented in GeomDrape tool [18]. This software provides a fibre quality chart (showing the fibre distortions, the rate of falling, and the rate of draped surface) to predict difficult impregnated regions.  It can be used to optimise the draping process (with respect to the above quality measure) by improving the layup directions or the marker data location.  The layup of complex curved surfaces can be made in a few seconds. The use of draping tool should not only allow for a more detailed assessment of the draped fabric (Figure 2(c)), but also for an optimization in terms of total fabric shear deformation (Figure 2(e)), flat pattern of the initial fabric (Figure 2(d)) and fibre stretching (Figure 2(f)), and the optimization of the best layup start point (Figure 2(a)).


## 2.2. Application of Geometrical Approach

Two draping simulation examples are given. These simulations are performed using the geometrical analysis computer code GeomDrap [18]. For each example, we assume that a mesh of the piece to drape is given. The first example shows the influence of the fibre orientations draping in the draping process. The second example shows the efficiency of the proposed method to simulate geometrically the draping composite part.The first example concerns the geometrical draping of a complex piece composed of two half hemispheres with a radius of 38.8 mm that are connected by a half cylinder with a length of 170 mm [23]. Figure 3 shows the piece as well as a mesh of the mechanical part. The centroid of the part is chosen as the impact point for which two different fibre orientations are specified. Figure 4 shows the resulting 3D surface layup for the 0°/90° fibre orientation (Figure 4(a)) and the corresponding 2D flat patterns (Figure 4(b)). Likewise, Figure 5 shows the draping results for the ±45° fibre orientation (Figure 5(a) 3D surface layup and Figure 5(b) the corresponding 2D flat pattern). The shears angles between fibres are given in Figure 6. One can notice that, in the considered cases, the surface of the piece is draped globally. However, to drape completely the proposed surface without excessive fibre distortion <50°, it is necessary to make cuts chisel along the direction 45° for 0°/90° fibre orientation (Figure 6(a)) and along the direction 0°  ± 45° fibre orientation (Figure 6(b)). In case we can allow excessive distortions of the fibers without cutting chisel, the results given in Figure 7 show the distribution of fibers in both directions 0°/90° (Figure 7(a)) and ±45° (Figure 7(a)). We can see that the shear limit reached 52° in the case of 0°/90° fibre orientation and 89° in the case of ±45° fibre orientation.Figure 3
CAD of complex part.Draping of (0°/90°) fabric: (a) 3D surface layup and (b) corresponding 2D flat pattern.
(a)
(b)Draping of (−45°/45°) fabric: (a) 3D surface layup and (b) corresponding 2D flat pattern.
(a)
(b)Shear angles: (a) 0°/90° fabric orientation and (b) −45°/+45° fabric orientation.
(a)
(b)Optimal 3D surface layup: (a) 0°/90° fabric orientation and (b) −45°/+45° fabric orientation.
(a)
(b)To compare the experimentally obtained fibre angle result [23] with the geometrical predicted ones, two cross sections where most shearing occurs are examined. The amount of shear along the edge of these cross sections is used to compare with the draping approach. Figure 3 gives an overview of the two cross sections along the symmetrical L1 and along the median line L2 for different ply orientations. The line L2 corresponds with the examined cross section for a 0°/90° fabric draping orientation and the line L1 for a 45 ± 45° fabric draping orientation. From Figure 8(a), it can be concluded that for the 0°/90° orientation, the agreement between the experimental shear angles and the geometrical draping results is good. On the other side, the drape orientation ±45°, the results from the geometrical model do not agree at all with the experimental results [23] (see Figure 8(b)). The oversimplification of the fabric deformation in this geometrical model gives shear angles up to 89° in case of the 45° fabric orientation, which is impossible for fabric composite. The geometrical model is used with a cutoff shear angle based either on an experimentally determined locking angle or the maximum orientation that the designer is prepared to tolerate, when defining a cutoff angle of 38°, which equals the experimentally determined locking angle. From the above discussed results, it can be concluded that the geometrical draping method is not appropriate for predicting the fiber directions of a fabric sheet after stamp forming.Shear angles of optimal geometrical draping: (a) along the line L2 of 0°/90° and (b) along the line L1 of −45°/45°.
(a)
(b)The second example is the draping of a complex shape (car hood). The centroid of the part is chosen as the starting point from which the 0°/90° and ±45° fibre orientations are specified. Figure9 shows the resulting 3D draping for the two orientations (Figure 9(a) for 0°/90° fabric orientation and Figure 9(b) for 45°/45° fabric orientation). We can note that all the part surface is completely draped. Figure 10 presents shaded contours interpolated from the map of the fiber distortions of 0°/90° (Figure 10(a)) and ±45° fiber orientations (Figure 10(b)). The fiber distortions for both 0°/90° and ±45° draping are very small but the maximum shear angle localization is different.3D geometrical draping of hood car part: (a) 0°/90° fabric orientation and (b) −45°/45° fabric orientation.
(a)
(0°/90°)
(b)
(−45°/45°)Isovalues of shear angles between fibres: (a) 0°/90° fabric orientation and (b) −45°/45° fabric orientation.
(a)
(0°/90°)
(b)
(−45°/45°)


## 3. Mechanical Finite Element Approach

The mechanical behaviour of composite fabric forming is complex, and it is a multiscale problem due to the fibres behavior and the interactions of the fibres or yarns [37–41]. The macroscopic behaviour is very much dependent on the interactions of fibres at the mesoscale (scale of the woven unit cell) and at the microscale (level of the fibres constituting yarns). Despite of the great amount of work in the field, there is no widely accepted model that accurately describes all the main aspects of fabric mechanical behaviour. The main model families come from the multiscale nature of the textile. During the forming process of woven fabric, the two main modes of deformation at the mesoscopic scale are the stretching of the fibres due fibres undulation and the in-plane shearing of the fabric resulting in a change of the angle between the warp and the left yarns. In the deep drawing or the draping of woven fabrics, the in-plane shear of fibres is the principal mode of deformation and is very different than the sheet metal [37, 40]. Figure 11 shows the evolution of two straight lines draw alternatively on warp and weft fibre directions during the forming deformation. These lines become curved but remain continuous. The absence of interyarn sliding ensured by the fabric weaving, viscoelastic behavior of resin, and friction fiber/fiber and fiber/resin can be observed over the main areas of the fabric (i.e., far enough from the free edges of the fabric). Also, for the composite fabrics based on high modulus, the compressive as well as bending stiffness is negligible compared to the in-plane membrane stiffness. The assumption is that each cross connection of straight warp and weft fibres before deformation remains cross connected during the deformation. The basic assumptions for the mechanical forming are that the woven fabric is considered as a continuous 3D material. The warp and weft fibres are represented by a truss to which connecting points are hinged, and the membrane resin is coupled kinematically to the fabric at these connecting points.Figure 11
Fibres deformation during the forming process.


### 3.1. Model Description

In the unit cell of the mesoscopic model used here for a plain weave fabric, each fibre and resin is modelled and is assumed to be a straight or a curved beam or truss. The tensile load is carried by the 3D truss or beam elements that will capture the changes in the orthotropic directions during the shearing. The 3D membrane or shell element accounts only for the shearing resistance of the fabric and hence has no tensile stiffness. The interaction warp-weft yarn and resin fibres is negligible. In this approach, the stress and strain of a continuous material are related to fibrous reinforcement using the constitutive relation in a nonorthogonal frame directed by the fibre directions. In this study, we consider two yarn directions, and we use them to define the nonorthogonal frame. This approach uses the Green-Naghdi frame. It is an orthonormal frame which is rotated byR, the rotation of the polar decomposition in which the local stress increment computations at finite strain are made [42–44].Two reference frames have to be considered. Theei unit vectors define the local orthogonal reference frame that rotates with the continuum material, and the gi basis vectors form a nonorthogonal frame that follows the fibre direction. Here, g1,g2 correspond to warp and weft directions, respectively. In this case, for each connecting point X→f of warp and weft yarns is associated a material position space of a resin X→m. At the connecting points, we have X→f=X→m=X→ before deformation. The current position of these points is obtained by
(1)dx→f=Ff(X→,t)dX→ffibres,dx→m=Fm(X→,t)dX→mresin,
where Ff and Fm are the deformation gradient tensor of fibre and resin, respectively. The relationship of the no sliding interfibre can be written at each connecting point as x→=x→(X→,t)/x→f=x→m=x→. The gradient of transformation of the preimpregnated woven fabric F and the pseudo gradient of transformation of the fibre Ff are defined by the function
(2)Fijf=λLgi⊗g0jfibres,Fijm=∂xm∂Xme0i⊗e0jresin,
where λLf is the longitudinal elongation of each fibre and g0i, and gi are, respectively, the fibre orientations in the initial C0 and the current Ct configurations and are the resin. The relative rotation of fibre can be associated to the rotation of the rigid body of the median line of the fibre R. We can note this assumption by the following kinematic relation:
(3)R=gi⊗g0i.
Using the previous assumptions, the mechanical deformation of composite fabric depends on the relative movement of fibres and the deformation of unpolymerized resin. Consider
(4)λif=g0iFfTFfg0iffibres,Um=FmTFmresin.
The shaping problem imposes the use of incremental formulation in finite deformations. In finite deformation analysis, a careful distinction has to be made between the coordinate systems that can be chosen to describe the behaviour of the body. The rate constitutive equations for finite strain use objective derivatives [43]. The problem of the integration of strain rate tensors is a central one in large deformations. The rate of deformation tensor of woven fabric is obtained by
(5)Df=(λ˙ifλif)(gi⊗gi)fibres,Dm=12(F˙F-1+F-TF˙T)(e0i⊗e0i)resin.
The rate equations for finite strains use objective derivatives [44]. The approaches traditionally developed in finite element codes for anisotropic metal at large strains are based on Jaumann corotational formulation or the Green-Naghdi approach. In these models, a rotation is used both to define an objective derivative for the hypoelastic model law and to update the orthotropic frame. The rotations used in Green-Naghdi and Jaumann derivatives are average rotations of the material (polar rotation and corotational rotations, respectively). The frame associated with Green-Naghdi’s derivative is defined, at the material point considered. The rotation Rf is used to update the initial constitutive axes of warps or weft fibres g0i⊗g0i to the current constitutive axes, and the rotation R is used to update the initial constitutive axes of resin e0i⊗e0i. The stretching tensors in (5) are written in the rigid body rotation frames. The longitudinal component is D-Lf=λ˙Lf/λLf, and the transversal components (D-Tf,D-3f) are obtained by the unidirectional behaviour of fibre as D-TfR=D-3fR=-νLTD-LfR⇒λTf=λ3f=(λLf)-νLT, and νLT is the Poison’s ratio of fibre. Using Green-Naghdi’s objective tensor stress, the stress rate σ-Lf, depending on the stretching deformation D-Lf, and the stress rate tensor of the membrane resin σ-m, depending on the tensor deformation rate D-m and the elastic properties Cm,  can be written at each time as follows. The stress and strain relationship of the woven fabric composites are defined in the nonorthogonal material coordinate frame. Therefore, coordinate transformations of stress and strain into the orthogonal coordinate system should be considered in the nonorthogonal constitutive model as σ=σ-f+σ-m as follows:(6)σ=(∑weftELf(λ1f)λ˙1fλ1f-νLTELf(λ2f)∑weftλ˙2fλ2f000∑warpELf(λ2f)λ˙2fλ2f-νLT∑weftELf(λ1f)λ˙1fλ1f0000)+(Em1-νm2D-11mνmEm1-νm2D-22m0νmEm1-νm2D-11mEm1-νm2D-22m000GmD-12m).The constitutive law of fibres is nonlinear and is written in terms of longitudinal modulus of stretching ELf(λLf), and the compressive stiffness (λLf≤0) of fiber is supposed negligible ELf=0. Later is the function of elongation of warp and weft fibre (λ1f,λ2f), effective elastic modulus of fibre E-f, and undulation factor εsh; Em,  νm are the membrane elastic properties. To determine the stress state in a membrane material at a given time, the deformation history must be considered. For linear viscoelastic materials, a superposition of hereditary integrals describes the time-dependent response. Let Gm(t) be the shear stress relaxation modulus of the unpolymerized resin and Gm∞=Gm(t=∞) the limit value. The viscoelastic behaviour of unpolymerized resin is formulated in the time domain by the hereditary integral and using the relaxation time τk and the shear modulus relaxation, which are material parameters Gmk. Hereditary integrals with Prony series kernels can be applied to model the shear behavior of the unpolymerized resin. The behavior of fibre and resin can be written as
(7)ELf(λLf)=E-f(1-exp(-λ˙LfλLfεsh))fibres,Gm(t)=Gm(∞)+∑1kGmkexp(-tτk)resin.


### 3.2. Finite Element Formulation

Each material point is moving as in a continuum, ensured by the nonsliding of fibers due to fabric weaving, and resin behavior. Therefore, a nodal approximation for the displacement can be used. The deformation of composite fabric is described within the frame of membrane assumptions. The energy of deformationΠ(u˙) is obtained by a summation of membrane strain energy of  unpolymerized resin and elastic tensile strain energy of fibers
(8)δΠ(u˙)=h0∫S0mσ-m:δD-mds+∑fibresS0f∫Lfσ-Lf:δD-Lfdl-∫Γσt-·δu˙dΓ,
where h0 denote the initial thickness of fabric, Lf the length of fiber, and S0f the initial effective cross section of fibre, and t- is the external surface load applied along Γσ of woven fabric. The effective cross section of the fibre S0f and the effective surface S0m of the membrane resin, that assumed no void between the fibre and the resin, were used and are calculated by using the fibre volume fraction Vf of the woven fabric (Vfabric=VfS0fLf+(1-Vf)h0S0m). The global equilibrium of the fabric is obtained by minimizing the total potential energy  Π(u˙). The effect of spatial equilibrium of composite material on the actual configuration is established in terms of nonlinear equations: kinematic nonlinearity, material nonlinearity, and contact with friction nonlinearity. It is linearized for each load increment by an iterative Newton method. It should be emphasized that during the motion, nodes and elements are permanently attached to the material points with which they were initially associated. Consequently, the subsequent motion is fully described in terms of the current nodal positions as
(9)x=∑k=1mNk(ξ,η)Xk+∑k=1mNk(ξ,η)uk,
where uk are the nodal displacements of each connecting point, Nk(ξ,η) are the standard shape functions (of membrane of truss element), and (m) denotes the number of nodes.The discretization of rate deformation tensor can be obtained by introducing (5) into the definition of D-Lf and D-m  and given in (8) to give
(10)D-Lf=λ˙LfλLf=1αn2(dxTdξdu˙dξ)=[Bfibres]u˙nfibres,D-m=12RT(du˙dξdξdx+∂u˙Tdξdξdx)R=[Bresin]u˙nresin,
where
(11)αn=dXfTdξdXfdξ+2dXfTdξdudξ+duTdξdudξ,
and [Bresine] and [Bfibrese] are the geometric or strain-displacement matrix of membrane resin and truss fibres.In a finite element approximation, the only independent variables in the equations of linearized virtual work are the displacements of the material points. Substituting the element coordinate and displacement interpolations into the equilibrium equations, for a given set of elements, we obtain that the state that forces acting on a fabric equals the mass times the acceleration of the body(12)∑e[Me]{u¨}+∑e({ℜinte}-{ℜexte})={0},
where [Me] is the consistent composite mass matrix, and ({ℜinte}-{ℜexte}) is the so-called quasistatic equilibrium residual
(13)[Me]=ρrh0∫Sm[Ne]T[Ne]ds+∑fibresρfS0f∫L[N~e]T[N~e]dl,{ℜinte}-{ℜexte}=h0∫Sm[Bresine]T{σ-mR}ds-∑fibresS0f∫LfE-Lf[Bfibrese]T{σ-fR}dl-∫Γf[Ne]T{t-}ds.[Ne] and [N~e] are the matrix of the nodal interpolation functions of both associated elements, and h0,S0f are the initial thickness and surface of resin and fibres, respectively. The index e refers to the eth element.According to the different modes of deformation occurring in the prepreg fabric during the shaping process, bicomponent finite elements are developed to characterize the mechanical behavior of thin composite structures. The bicomponent element is based on an association of 3D linear membrane finite elements (T3 and Q4) combined with a complementary truss linear finite elements. The global stiffness of composite fabric is obtained by the summation of elementary stiffness matrix of warp fiber, elementary stiffness matrix of weft fiber, and elementary stiffness matrix of resin. These finite elements are complementary in the finite element discretization, isoparametric and use three dof per node, and use the same mechanical formulation in finite deformations (Green-Naghdi’s approach). The nonlinear constitutive equation of fibre behavior is implemented in the Abaqus/Explicit using VUMAT user’s subroutine [13–15]. The governing equilibrium equation (12) is solved as a dynamic problem using explicit integration. This is achieved by using the central difference method to approximate the velocity and the acceleration in the next time step, using only information from the previous step, where all state variables are known. This approach has proven to be, in particular, suitable to highly nonlinear geometric and material problems, particularly where a large amount of contact between different structural parts occurs. In the present work, the dynamic explicit (DE) resolution procedure is used within the general purpose FE code Abaqus/Explicit. The DE algorithm available in Abaqus/Explicit for solving the algebraic system works by using the lumped form of the mass matrix [26]. The major disadvantage of the explicit scheme is that it is only stable for time steps small enough. This can cause, for instance, the energy balance to be changed. One simple and conservative criterion for a stable time increment is
(14)tstable=ρrLeE-f,
where Le is the minimal length of membrane element, ρr is the resin density, and E-f is Young’s modulus of fibre. The interpretation of this requirement is that the time increment must be shorter than the time it takes for a propagating wave in the material to cross the shortest side of an element. Thus, the element size and the critical time increment are connected. The smaller the elements are, the shorter the time increment must be. Because of the difference in wave speed for different materials, the critical element size is larger for a stiff material (e.g., steel) than for a softer material (e.g., polymer).


### 3.3. Application of Mechanical Finite Element

#### 3.3.1. Uniaxial Test

Due to the importance of the composite fabric behaviour on material formability, tensile test of preimpregnated fabric is proposed in order to study the influence of fibre orientation, fibre undulations due to fabric weaving and resin behaviour. At processing, the impregnated fabric is idealized as a viscous material subject to the kinematic constraints of incompressibility and inextensibility in the fibre direction. The preimpregnated fabric tested in this study was a satin 5 with aramid woven fabric. The fabric was impregnated with epoxy resin using a hot-melt prepregging process (the mechanical properties are given in Table1).  A lower loading velocity will generate lower viscous forces at the intraply shearing of fabric. In this experiment, a displacement in the vertical direction is imposed at the moving extremity of the rectangular composite specimen with three layers of fibres (length = 150 mm, width = 30 mm, and thickness = 2 mm). The uniaxial tensile test (see Figure 12) is carried out for different orientations of fibres with the loading direction (0, 15, 30, and 45°). The experimental effort imposed by the tensile machine is compared to the numerical values for different fibre orientations in Figure 13. In this figure, we can show the good correlation between the model and the experimental results. The agreement between predicted and experimental values is good and proves the validity of the proposed model of preimpregnated woven fabric behaviour. The numerical model described previously clearly shows the strong nonlinearity of this behaviour law. It takes into account the mechanical characteristics of a viscoelastic resin, the anisotropic behaviour of fabric, and the geometrical nonlinearities due to the high deformability of fibres (straightening and relative rotation).Table 1
Mechanical properties of the pre-impregnated composite fabric.
(a)
E
f (MPa)
ε
sh
ρ (g/cm3)
130000
0.005
1,45
(b)
Time (s)
Shear modulus of resinGmk
0.01
0.02332
0.1
0.023332
1
0.083509
10
0.11723
100
0.14423
1000
0.178Tensile test of preimpregnated woven fabric.
(a)
(b)Figure 13
Load force versus displacement for different fibre orientations.


#### 3.3.2. Deep Drawing of Dry Woven Fabric

The second example is the 3D deep drawing of dry glass woven fabric by hemispherical punch by experimental and mechanical approaches. The initial shape of the glass fibre fabric is a square (360×360 mm). Its edges are free, but a pressure equal to 2 MPa is applied on the binder and the friction coefficient between the glass fabric and the steel tools is 0.27. The forming simulation has been performed in Abaqus/Explicit using rigid tools. As mentioned before, the shear resistance model is assigned to reduced membrane elements, and truss elements parallel to the membrane edges represent the high tensile stiffness in the yarn direction. Figures 14 and 15 report, respectively, the experimental obtained shapes with respect to 0°/90° (Figure 14(a)) and ±45° fibre orientations (Figure 14(a)) and the corresponding final shape using the numerical approach (Figure 15(a) for 0°/90° and Figure 15(b)  for ±45° fibre orientations). Likewise, the numerical simulation agrees with the experimental results. The angular fibre distortion exceeds 38° along the diagonal axis for 0°/90° fibre orientations and 50° along the diagonal axis for ±45° fibre orientations. Furthermore, we can notice that the final shape obtained with 0°/90° fibre direction is very different from the −45°/+45°. Another interesting result of the numerical calculation is the angular distortion variations between warp and weft fibres.Experimental shape: (a) 0°/90° fibre orientations and (b) ±45° fibre orientations.
(a)
(b)Predicted shape and shear angles: (a) 0°/90° and (b) ±45° fibre orientations.
(a)
(b)


#### 3.3.3. Deep Drawing of Prepreg Woven Fabric

The numerical analysis of prepreg composite fabric deformation by the deep drawing process is performed by utilizing the commercial FEM package Abaqus/Explicit. The resin is modeled by using 1600 membrane finite elements (linear triangular element M3D3), and warp and weft fibres are modeled by 3200 truss finite elements (linear element T3D2). The rigid surface is modeled by 1600 Bezier patches (three nodes R3D3 and four nodes R3D4). The behavior of the resin is assumed to be isotropic viscoelastic, and the behaviour of the fibre is supposed to be elastic.The first forming example concerns the 3D deep drawing of aramid preimpregnated fabric with conical tools (see Figure16). The mechanical properties of the used material are given in  Table 1.  Figure 17 reports the experimentally obtained shapes with respect to 0°/90° (Figure 17(a)) and ±45° fibre orientation (Figure 17(b)) for different punch displacements (initial, and final for 100 mm of punch displacement). The corresponding final computed shapes are shown in Figure 18(a) for 0°/90° and in Figure 18(b) for ±45° fibre orientations. The corresponding final geometrical shapes are shown in Figure 18(c) for 0°/90° and in Figure 18(d) for ±45° fibre orientations. There is good agreement between the experimental and the predicted (mechanical and geometrical) results.Figure 16
Geometry of deep drawing tools.Experimental results: (a) 0°/90° fibre orientation and (b) −45°/45° fibre orientation.
(a)
(0°/90°)
(b)
(−45°/+45°)Geometrical and FE predicted results 0°/90° and −45°/45° fibre orientations.
(a)
FEM (0°/90°)
(b)
FEM (−45°/45°)
(c)
GEOMETRICAL  (0°/90°)
(d)
GEOMETRICAL (−45°/45°)We notice that these shear angle values are very large >38° for mechanical approach along the median line for ±45° fabric (Figure19(b)) and along the diagonal line for 0°/90° fabric (Figure 19(a)). But along the median lines of 0°/90° and along the diagonal lines of ±45°, the angular distortions are very small <6°. There is good agreement between the geometrical and the mechanical modeling.Shear angles: (a) along the diagonal line of 0°/90° and (b) along the symmetrical line of 45°/45°.
(a)
(b)


#### 3.3.4. Deep Drawing of Composite with Complex Punch

The third forming concerns the simulation of composite glass fabric stamping with complex tools using the mechanical approach and the computational remeshing procedure. The initial shape of the taffetas woven fabric is700×350 mm. The high tensile stiffness along the warp and weft yarn direction, respectively, is introduced via truss elements that connect the nodes of the membrane element. The Young’s modulus in warp and weft directions is 70 GPa. The local fibre directions are tracked on the part via a 3D-digital image correlation technique that assesses the intersections of a raster pattern that is painted on the preconsolidated part along the fibres directions. The geometry of the stamping simulation is shown in Figure 20. The first order triangular T3 and quadrilateral Q4 membrane element types were used for the stamping. The binding force was 30 N, and the punch stroke was 47 mm. The friction coefficient between the glass fabric and the tools is 0.3. The amounts of draw-in and fiber orientation were measured. The fibre direction is assumed to be 0°/90°. The final deformed shape of 0°/90° fabric for a 47 mm displacement of the punch is shown in Figure 21(a) and is compared to the experimental one (Figure 21(b)). The shear angle distribution was plotted in Figure 22. We can note the adapting meshes of the warp and the weft to the curved shape of the punch. As shown in the Figure 23, maximum shear angle was observed at the highly curved area, and such positions were almost the same for all simulation cases. From the numerical simulation, maximum shear angle was 35°.CAD and experimental tools of composite stamping.(a) Deformed mesh of the complex part and (b) experimental result.
(a)
(b)Figure 22
Isovalues of shear angles.Figure 23
Experimental and predicted shear angles along the diagonal line AB.


## 3.1. Model Description

In the unit cell of the mesoscopic model used here for a plain weave fabric, each fibre and resin is modelled and is assumed to be a straight or a curved beam or truss. The tensile load is carried by the 3D truss or beam elements that will capture the changes in the orthotropic directions during the shearing. The 3D membrane or shell element accounts only for the shearing resistance of the fabric and hence has no tensile stiffness. The interaction warp-weft yarn and resin fibres is negligible. In this approach, the stress and strain of a continuous material are related to fibrous reinforcement using the constitutive relation in a nonorthogonal frame directed by the fibre directions. In this study, we consider two yarn directions, and we use them to define the nonorthogonal frame. This approach uses the Green-Naghdi frame. It is an orthonormal frame which is rotated byR, the rotation of the polar decomposition in which the local stress increment computations at finite strain are made [42–44].Two reference frames have to be considered. Theei unit vectors define the local orthogonal reference frame that rotates with the continuum material, and the gi basis vectors form a nonorthogonal frame that follows the fibre direction. Here, g1,g2 correspond to warp and weft directions, respectively. In this case, for each connecting point X→f of warp and weft yarns is associated a material position space of a resin X→m. At the connecting points, we have X→f=X→m=X→ before deformation. The current position of these points is obtained by
(1)dx→f=Ff(X→,t)dX→ffibres,dx→m=Fm(X→,t)dX→mresin,
where Ff and Fm are the deformation gradient tensor of fibre and resin, respectively. The relationship of the no sliding interfibre can be written at each connecting point as x→=x→(X→,t)/x→f=x→m=x→. The gradient of transformation of the preimpregnated woven fabric F and the pseudo gradient of transformation of the fibre Ff are defined by the function
(2)Fijf=λLgi⊗g0jfibres,Fijm=∂xm∂Xme0i⊗e0jresin,
where λLf is the longitudinal elongation of each fibre and g0i, and gi are, respectively, the fibre orientations in the initial C0 and the current Ct configurations and are the resin. The relative rotation of fibre can be associated to the rotation of the rigid body of the median line of the fibre R. We can note this assumption by the following kinematic relation:
(3)R=gi⊗g0i.
Using the previous assumptions, the mechanical deformation of composite fabric depends on the relative movement of fibres and the deformation of unpolymerized resin. Consider
(4)λif=g0iFfTFfg0iffibres,Um=FmTFmresin.
The shaping problem imposes the use of incremental formulation in finite deformations. In finite deformation analysis, a careful distinction has to be made between the coordinate systems that can be chosen to describe the behaviour of the body. The rate constitutive equations for finite strain use objective derivatives [43]. The problem of the integration of strain rate tensors is a central one in large deformations. The rate of deformation tensor of woven fabric is obtained by
(5)Df=(λ˙ifλif)(gi⊗gi)fibres,Dm=12(F˙F-1+F-TF˙T)(e0i⊗e0i)resin.
The rate equations for finite strains use objective derivatives [44]. The approaches traditionally developed in finite element codes for anisotropic metal at large strains are based on Jaumann corotational formulation or the Green-Naghdi approach. In these models, a rotation is used both to define an objective derivative for the hypoelastic model law and to update the orthotropic frame. The rotations used in Green-Naghdi and Jaumann derivatives are average rotations of the material (polar rotation and corotational rotations, respectively). The frame associated with Green-Naghdi’s derivative is defined, at the material point considered. The rotation Rf is used to update the initial constitutive axes of warps or weft fibres g0i⊗g0i to the current constitutive axes, and the rotation R is used to update the initial constitutive axes of resin e0i⊗e0i. The stretching tensors in (5) are written in the rigid body rotation frames. The longitudinal component is D-Lf=λ˙Lf/λLf, and the transversal components (D-Tf,D-3f) are obtained by the unidirectional behaviour of fibre as D-TfR=D-3fR=-νLTD-LfR⇒λTf=λ3f=(λLf)-νLT, and νLT is the Poison’s ratio of fibre. Using Green-Naghdi’s objective tensor stress, the stress rate σ-Lf, depending on the stretching deformation D-Lf, and the stress rate tensor of the membrane resin σ-m, depending on the tensor deformation rate D-m and the elastic properties Cm,  can be written at each time as follows. The stress and strain relationship of the woven fabric composites are defined in the nonorthogonal material coordinate frame. Therefore, coordinate transformations of stress and strain into the orthogonal coordinate system should be considered in the nonorthogonal constitutive model as σ=σ-f+σ-m as follows:(6)σ=(∑weftELf(λ1f)λ˙1fλ1f-νLTELf(λ2f)∑weftλ˙2fλ2f000∑warpELf(λ2f)λ˙2fλ2f-νLT∑weftELf(λ1f)λ˙1fλ1f0000)+(Em1-νm2D-11mνmEm1-νm2D-22m0νmEm1-νm2D-11mEm1-νm2D-22m000GmD-12m).The constitutive law of fibres is nonlinear and is written in terms of longitudinal modulus of stretching ELf(λLf), and the compressive stiffness (λLf≤0) of fiber is supposed negligible ELf=0. Later is the function of elongation of warp and weft fibre (λ1f,λ2f), effective elastic modulus of fibre E-f, and undulation factor εsh; Em,  νm are the membrane elastic properties. To determine the stress state in a membrane material at a given time, the deformation history must be considered. For linear viscoelastic materials, a superposition of hereditary integrals describes the time-dependent response. Let Gm(t) be the shear stress relaxation modulus of the unpolymerized resin and Gm∞=Gm(t=∞) the limit value. The viscoelastic behaviour of unpolymerized resin is formulated in the time domain by the hereditary integral and using the relaxation time τk and the shear modulus relaxation, which are material parameters Gmk. Hereditary integrals with Prony series kernels can be applied to model the shear behavior of the unpolymerized resin. The behavior of fibre and resin can be written as
(7)ELf(λLf)=E-f(1-exp(-λ˙LfλLfεsh))fibres,Gm(t)=Gm(∞)+∑1kGmkexp(-tτk)resin.


## 3.2. Finite Element Formulation

Each material point is moving as in a continuum, ensured by the nonsliding of fibers due to fabric weaving, and resin behavior. Therefore, a nodal approximation for the displacement can be used. The deformation of composite fabric is described within the frame of membrane assumptions. The energy of deformationΠ(u˙) is obtained by a summation of membrane strain energy of  unpolymerized resin and elastic tensile strain energy of fibers
(8)δΠ(u˙)=h0∫S0mσ-m:δD-mds+∑fibresS0f∫Lfσ-Lf:δD-Lfdl-∫Γσt-·δu˙dΓ,
where h0 denote the initial thickness of fabric, Lf the length of fiber, and S0f the initial effective cross section of fibre, and t- is the external surface load applied along Γσ of woven fabric. The effective cross section of the fibre S0f and the effective surface S0m of the membrane resin, that assumed no void between the fibre and the resin, were used and are calculated by using the fibre volume fraction Vf of the woven fabric (Vfabric=VfS0fLf+(1-Vf)h0S0m). The global equilibrium of the fabric is obtained by minimizing the total potential energy  Π(u˙). The effect of spatial equilibrium of composite material on the actual configuration is established in terms of nonlinear equations: kinematic nonlinearity, material nonlinearity, and contact with friction nonlinearity. It is linearized for each load increment by an iterative Newton method. It should be emphasized that during the motion, nodes and elements are permanently attached to the material points with which they were initially associated. Consequently, the subsequent motion is fully described in terms of the current nodal positions as
(9)x=∑k=1mNk(ξ,η)Xk+∑k=1mNk(ξ,η)uk,
where uk are the nodal displacements of each connecting point, Nk(ξ,η) are the standard shape functions (of membrane of truss element), and (m) denotes the number of nodes.The discretization of rate deformation tensor can be obtained by introducing (5) into the definition of D-Lf and D-m  and given in (8) to give
(10)D-Lf=λ˙LfλLf=1αn2(dxTdξdu˙dξ)=[Bfibres]u˙nfibres,D-m=12RT(du˙dξdξdx+∂u˙Tdξdξdx)R=[Bresin]u˙nresin,
where
(11)αn=dXfTdξdXfdξ+2dXfTdξdudξ+duTdξdudξ,
and [Bresine] and [Bfibrese] are the geometric or strain-displacement matrix of membrane resin and truss fibres.In a finite element approximation, the only independent variables in the equations of linearized virtual work are the displacements of the material points. Substituting the element coordinate and displacement interpolations into the equilibrium equations, for a given set of elements, we obtain that the state that forces acting on a fabric equals the mass times the acceleration of the body(12)∑e[Me]{u¨}+∑e({ℜinte}-{ℜexte})={0},
where [Me] is the consistent composite mass matrix, and ({ℜinte}-{ℜexte}) is the so-called quasistatic equilibrium residual
(13)[Me]=ρrh0∫Sm[Ne]T[Ne]ds+∑fibresρfS0f∫L[N~e]T[N~e]dl,{ℜinte}-{ℜexte}=h0∫Sm[Bresine]T{σ-mR}ds-∑fibresS0f∫LfE-Lf[Bfibrese]T{σ-fR}dl-∫Γf[Ne]T{t-}ds.[Ne] and [N~e] are the matrix of the nodal interpolation functions of both associated elements, and h0,S0f are the initial thickness and surface of resin and fibres, respectively. The index e refers to the eth element.According to the different modes of deformation occurring in the prepreg fabric during the shaping process, bicomponent finite elements are developed to characterize the mechanical behavior of thin composite structures. The bicomponent element is based on an association of 3D linear membrane finite elements (T3 and Q4) combined with a complementary truss linear finite elements. The global stiffness of composite fabric is obtained by the summation of elementary stiffness matrix of warp fiber, elementary stiffness matrix of weft fiber, and elementary stiffness matrix of resin. These finite elements are complementary in the finite element discretization, isoparametric and use three dof per node, and use the same mechanical formulation in finite deformations (Green-Naghdi’s approach). The nonlinear constitutive equation of fibre behavior is implemented in the Abaqus/Explicit using VUMAT user’s subroutine [13–15]. The governing equilibrium equation (12) is solved as a dynamic problem using explicit integration. This is achieved by using the central difference method to approximate the velocity and the acceleration in the next time step, using only information from the previous step, where all state variables are known. This approach has proven to be, in particular, suitable to highly nonlinear geometric and material problems, particularly where a large amount of contact between different structural parts occurs. In the present work, the dynamic explicit (DE) resolution procedure is used within the general purpose FE code Abaqus/Explicit. The DE algorithm available in Abaqus/Explicit for solving the algebraic system works by using the lumped form of the mass matrix [26]. The major disadvantage of the explicit scheme is that it is only stable for time steps small enough. This can cause, for instance, the energy balance to be changed. One simple and conservative criterion for a stable time increment is
(14)tstable=ρrLeE-f,
where Le is the minimal length of membrane element, ρr is the resin density, and E-f is Young’s modulus of fibre. The interpretation of this requirement is that the time increment must be shorter than the time it takes for a propagating wave in the material to cross the shortest side of an element. Thus, the element size and the critical time increment are connected. The smaller the elements are, the shorter the time increment must be. Because of the difference in wave speed for different materials, the critical element size is larger for a stiff material (e.g., steel) than for a softer material (e.g., polymer).


## 3.3. Application of Mechanical Finite Element

### 3.3.1. Uniaxial Test

Due to the importance of the composite fabric behaviour on material formability, tensile test of preimpregnated fabric is proposed in order to study the influence of fibre orientation, fibre undulations due to fabric weaving and resin behaviour. At processing, the impregnated fabric is idealized as a viscous material subject to the kinematic constraints of incompressibility and inextensibility in the fibre direction. The preimpregnated fabric tested in this study was a satin 5 with aramid woven fabric. The fabric was impregnated with epoxy resin using a hot-melt prepregging process (the mechanical properties are given in Table1).  A lower loading velocity will generate lower viscous forces at the intraply shearing of fabric. In this experiment, a displacement in the vertical direction is imposed at the moving extremity of the rectangular composite specimen with three layers of fibres (length = 150 mm, width = 30 mm, and thickness = 2 mm). The uniaxial tensile test (see Figure 12) is carried out for different orientations of fibres with the loading direction (0, 15, 30, and 45°). The experimental effort imposed by the tensile machine is compared to the numerical values for different fibre orientations in Figure 13. In this figure, we can show the good correlation between the model and the experimental results. The agreement between predicted and experimental values is good and proves the validity of the proposed model of preimpregnated woven fabric behaviour. The numerical model described previously clearly shows the strong nonlinearity of this behaviour law. It takes into account the mechanical characteristics of a viscoelastic resin, the anisotropic behaviour of fabric, and the geometrical nonlinearities due to the high deformability of fibres (straightening and relative rotation).Table 1
Mechanical properties of the pre-impregnated composite fabric.
(a)
E
f (MPa)
ε
sh
ρ (g/cm3)
130000
0.005
1,45
(b)
Time (s)
Shear modulus of resinGmk
0.01
0.02332
0.1
0.023332
1
0.083509
10
0.11723
100
0.14423
1000
0.178Tensile test of preimpregnated woven fabric.
(a)
(b)Figure 13
Load force versus displacement for different fibre orientations.


### 3.3.2. Deep Drawing of Dry Woven Fabric

The second example is the 3D deep drawing of dry glass woven fabric by hemispherical punch by experimental and mechanical approaches. The initial shape of the glass fibre fabric is a square (360×360 mm). Its edges are free, but a pressure equal to 2 MPa is applied on the binder and the friction coefficient between the glass fabric and the steel tools is 0.27. The forming simulation has been performed in Abaqus/Explicit using rigid tools. As mentioned before, the shear resistance model is assigned to reduced membrane elements, and truss elements parallel to the membrane edges represent the high tensile stiffness in the yarn direction. Figures 14 and 15 report, respectively, the experimental obtained shapes with respect to 0°/90° (Figure 14(a)) and ±45° fibre orientations (Figure 14(a)) and the corresponding final shape using the numerical approach (Figure 15(a) for 0°/90° and Figure 15(b)  for ±45° fibre orientations). Likewise, the numerical simulation agrees with the experimental results. The angular fibre distortion exceeds 38° along the diagonal axis for 0°/90° fibre orientations and 50° along the diagonal axis for ±45° fibre orientations. Furthermore, we can notice that the final shape obtained with 0°/90° fibre direction is very different from the −45°/+45°. Another interesting result of the numerical calculation is the angular distortion variations between warp and weft fibres.Experimental shape: (a) 0°/90° fibre orientations and (b) ±45° fibre orientations.
(a)
(b)Predicted shape and shear angles: (a) 0°/90° and (b) ±45° fibre orientations.
(a)
(b)


### 3.3.3. Deep Drawing of Prepreg Woven Fabric

The numerical analysis of prepreg composite fabric deformation by the deep drawing process is performed by utilizing the commercial FEM package Abaqus/Explicit. The resin is modeled by using 1600 membrane finite elements (linear triangular element M3D3), and warp and weft fibres are modeled by 3200 truss finite elements (linear element T3D2). The rigid surface is modeled by 1600 Bezier patches (three nodes R3D3 and four nodes R3D4). The behavior of the resin is assumed to be isotropic viscoelastic, and the behaviour of the fibre is supposed to be elastic.The first forming example concerns the 3D deep drawing of aramid preimpregnated fabric with conical tools (see Figure16). The mechanical properties of the used material are given in  Table 1.  Figure 17 reports the experimentally obtained shapes with respect to 0°/90° (Figure 17(a)) and ±45° fibre orientation (Figure 17(b)) for different punch displacements (initial, and final for 100 mm of punch displacement). The corresponding final computed shapes are shown in Figure 18(a) for 0°/90° and in Figure 18(b) for ±45° fibre orientations. The corresponding final geometrical shapes are shown in Figure 18(c) for 0°/90° and in Figure 18(d) for ±45° fibre orientations. There is good agreement between the experimental and the predicted (mechanical and geometrical) results.Figure 16
Geometry of deep drawing tools.Experimental results: (a) 0°/90° fibre orientation and (b) −45°/45° fibre orientation.
(a)
(0°/90°)
(b)
(−45°/+45°)Geometrical and FE predicted results 0°/90° and −45°/45° fibre orientations.
(a)
FEM (0°/90°)
(b)
FEM (−45°/45°)
(c)
GEOMETRICAL  (0°/90°)
(d)
GEOMETRICAL (−45°/45°)We notice that these shear angle values are very large >38° for mechanical approach along the median line for ±45° fabric (Figure19(b)) and along the diagonal line for 0°/90° fabric (Figure 19(a)). But along the median lines of 0°/90° and along the diagonal lines of ±45°, the angular distortions are very small <6°. There is good agreement between the geometrical and the mechanical modeling.Shear angles: (a) along the diagonal line of 0°/90° and (b) along the symmetrical line of 45°/45°.
(a)
(b)


### 3.3.4. Deep Drawing of Composite with Complex Punch

The third forming concerns the simulation of composite glass fabric stamping with complex tools using the mechanical approach and the computational remeshing procedure. The initial shape of the taffetas woven fabric is700×350 mm. The high tensile stiffness along the warp and weft yarn direction, respectively, is introduced via truss elements that connect the nodes of the membrane element. The Young’s modulus in warp and weft directions is 70 GPa. The local fibre directions are tracked on the part via a 3D-digital image correlation technique that assesses the intersections of a raster pattern that is painted on the preconsolidated part along the fibres directions. The geometry of the stamping simulation is shown in Figure 20. The first order triangular T3 and quadrilateral Q4 membrane element types were used for the stamping. The binding force was 30 N, and the punch stroke was 47 mm. The friction coefficient between the glass fabric and the tools is 0.3. The amounts of draw-in and fiber orientation were measured. The fibre direction is assumed to be 0°/90°. The final deformed shape of 0°/90° fabric for a 47 mm displacement of the punch is shown in Figure 21(a) and is compared to the experimental one (Figure 21(b)). The shear angle distribution was plotted in Figure 22. We can note the adapting meshes of the warp and the weft to the curved shape of the punch. As shown in the Figure 23, maximum shear angle was observed at the highly curved area, and such positions were almost the same for all simulation cases. From the numerical simulation, maximum shear angle was 35°.CAD and experimental tools of composite stamping.(a) Deformed mesh of the complex part and (b) experimental result.
(a)
(b)Figure 22
Isovalues of shear angles.Figure 23
Experimental and predicted shear angles along the diagonal line AB.


## 3.3.1. Uniaxial Test

Due to the importance of the composite fabric behaviour on material formability, tensile test of preimpregnated fabric is proposed in order to study the influence of fibre orientation, fibre undulations due to fabric weaving and resin behaviour. At processing, the impregnated fabric is idealized as a viscous material subject to the kinematic constraints of incompressibility and inextensibility in the fibre direction. The preimpregnated fabric tested in this study was a satin 5 with aramid woven fabric. The fabric was impregnated with epoxy resin using a hot-melt prepregging process (the mechanical properties are given in Table1).  A lower loading velocity will generate lower viscous forces at the intraply shearing of fabric. In this experiment, a displacement in the vertical direction is imposed at the moving extremity of the rectangular composite specimen with three layers of fibres (length = 150 mm, width = 30 mm, and thickness = 2 mm). The uniaxial tensile test (see Figure 12) is carried out for different orientations of fibres with the loading direction (0, 15, 30, and 45°). The experimental effort imposed by the tensile machine is compared to the numerical values for different fibre orientations in Figure 13. In this figure, we can show the good correlation between the model and the experimental results. The agreement between predicted and experimental values is good and proves the validity of the proposed model of preimpregnated woven fabric behaviour. The numerical model described previously clearly shows the strong nonlinearity of this behaviour law. It takes into account the mechanical characteristics of a viscoelastic resin, the anisotropic behaviour of fabric, and the geometrical nonlinearities due to the high deformability of fibres (straightening and relative rotation).Table 1
Mechanical properties of the pre-impregnated composite fabric.
(a)
E
f (MPa)
ε
sh
ρ (g/cm3)
130000
0.005
1,45
(b)
Time (s)
Shear modulus of resinGmk
0.01
0.02332
0.1
0.023332
1
0.083509
10
0.11723
100
0.14423
1000
0.178Tensile test of preimpregnated woven fabric.
(a)
(b)Figure 13
Load force versus displacement for different fibre orientations.


## 3.3.2. Deep Drawing of Dry Woven Fabric

The second example is the 3D deep drawing of dry glass woven fabric by hemispherical punch by experimental and mechanical approaches. The initial shape of the glass fibre fabric is a square (360×360 mm). Its edges are free, but a pressure equal to 2 MPa is applied on the binder and the friction coefficient between the glass fabric and the steel tools is 0.27. The forming simulation has been performed in Abaqus/Explicit using rigid tools. As mentioned before, the shear resistance model is assigned to reduced membrane elements, and truss elements parallel to the membrane edges represent the high tensile stiffness in the yarn direction. Figures 14 and 15 report, respectively, the experimental obtained shapes with respect to 0°/90° (Figure 14(a)) and ±45° fibre orientations (Figure 14(a)) and the corresponding final shape using the numerical approach (Figure 15(a) for 0°/90° and Figure 15(b)  for ±45° fibre orientations). Likewise, the numerical simulation agrees with the experimental results. The angular fibre distortion exceeds 38° along the diagonal axis for 0°/90° fibre orientations and 50° along the diagonal axis for ±45° fibre orientations. Furthermore, we can notice that the final shape obtained with 0°/90° fibre direction is very different from the −45°/+45°. Another interesting result of the numerical calculation is the angular distortion variations between warp and weft fibres.Experimental shape: (a) 0°/90° fibre orientations and (b) ±45° fibre orientations.
(a)
(b)Predicted shape and shear angles: (a) 0°/90° and (b) ±45° fibre orientations.
(a)
(b)


## 3.3.3. Deep Drawing of Prepreg Woven Fabric

The numerical analysis of prepreg composite fabric deformation by the deep drawing process is performed by utilizing the commercial FEM package Abaqus/Explicit. The resin is modeled by using 1600 membrane finite elements (linear triangular element M3D3), and warp and weft fibres are modeled by 3200 truss finite elements (linear element T3D2). The rigid surface is modeled by 1600 Bezier patches (three nodes R3D3 and four nodes R3D4). The behavior of the resin is assumed to be isotropic viscoelastic, and the behaviour of the fibre is supposed to be elastic.The first forming example concerns the 3D deep drawing of aramid preimpregnated fabric with conical tools (see Figure16). The mechanical properties of the used material are given in  Table 1.  Figure 17 reports the experimentally obtained shapes with respect to 0°/90° (Figure 17(a)) and ±45° fibre orientation (Figure 17(b)) for different punch displacements (initial, and final for 100 mm of punch displacement). The corresponding final computed shapes are shown in Figure 18(a) for 0°/90° and in Figure 18(b) for ±45° fibre orientations. The corresponding final geometrical shapes are shown in Figure 18(c) for 0°/90° and in Figure 18(d) for ±45° fibre orientations. There is good agreement between the experimental and the predicted (mechanical and geometrical) results.Figure 16
Geometry of deep drawing tools.Experimental results: (a) 0°/90° fibre orientation and (b) −45°/45° fibre orientation.
(a)
(0°/90°)
(b)
(−45°/+45°)Geometrical and FE predicted results 0°/90° and −45°/45° fibre orientations.
(a)
FEM (0°/90°)
(b)
FEM (−45°/45°)
(c)
GEOMETRICAL  (0°/90°)
(d)
GEOMETRICAL (−45°/45°)We notice that these shear angle values are very large >38° for mechanical approach along the median line for ±45° fabric (Figure19(b)) and along the diagonal line for 0°/90° fabric (Figure 19(a)). But along the median lines of 0°/90° and along the diagonal lines of ±45°, the angular distortions are very small <6°. There is good agreement between the geometrical and the mechanical modeling.Shear angles: (a) along the diagonal line of 0°/90° and (b) along the symmetrical line of 45°/45°.
(a)
(b)


## 3.3.4. Deep Drawing of Composite with Complex Punch

The third forming concerns the simulation of composite glass fabric stamping with complex tools using the mechanical approach and the computational remeshing procedure. The initial shape of the taffetas woven fabric is700×350 mm. The high tensile stiffness along the warp and weft yarn direction, respectively, is introduced via truss elements that connect the nodes of the membrane element. The Young’s modulus in warp and weft directions is 70 GPa. The local fibre directions are tracked on the part via a 3D-digital image correlation technique that assesses the intersections of a raster pattern that is painted on the preconsolidated part along the fibres directions. The geometry of the stamping simulation is shown in Figure 20. The first order triangular T3 and quadrilateral Q4 membrane element types were used for the stamping. The binding force was 30 N, and the punch stroke was 47 mm. The friction coefficient between the glass fabric and the tools is 0.3. The amounts of draw-in and fiber orientation were measured. The fibre direction is assumed to be 0°/90°. The final deformed shape of 0°/90° fabric for a 47 mm displacement of the punch is shown in Figure 21(a) and is compared to the experimental one (Figure 21(b)). The shear angle distribution was plotted in Figure 22. We can note the adapting meshes of the warp and the weft to the curved shape of the punch. As shown in the Figure 23, maximum shear angle was observed at the highly curved area, and such positions were almost the same for all simulation cases. From the numerical simulation, maximum shear angle was 35°.CAD and experimental tools of composite stamping.(a) Deformed mesh of the complex part and (b) experimental result.
(a)
(b)Figure 22
Isovalues of shear angles.Figure 23
Experimental and predicted shear angles along the diagonal line AB.


## 4. Conclusion

This paper has investigated the development of an efficient finite element and geometrical models for simulating the deformation of woven fabric preform materials under the influence of stitching patterns. The objective of this work is to show the influence of the woven fabric orientation on the prediction of the final shape of the woven reinforcement after forming using mechanical and geometrical approaches. The mechanical approach is based on a mesostructural model. It allows us to take into account the mechanical properties of fibres and resin and the various dominating modes of deformation of woven fabrics during the forming process. The geometrical approach is based on a modified Mosaic algorithm, which is suitable to generate a regular quad mesh representing the layup of the curved surfaces. Numerical examples concerning the deep drawing or draping of woven reinforcement demonstrated the efficiency of the proposed models.



---

*Source: 709495-2013-01-03.xml*