N Naouar

P Boisse

Kinematic modeling of transverse shear in textile composite reinforcements forming

  L'archive ouverte pluridisciplinaire

Introduction

Due to their outstanding performance (e.g., high specific stiffness, strength, and high designability), fiber-reinforced composites have been increasingly used for structural components in different fields such as transportation and civil engineering. The Liquid Composite Molding (LCM) process [START_REF] Sozer | Resin transfer molding (RTM) in polymer matrix composites[END_REF][START_REF] Deléglise | Modeling of high speed RTM injection with highly reactive resin with on-line mixing[END_REF] provides a possible option for the fabrication of composite components with high volume. The first step of the LCM process is forming dry fabric into the desired shape in molds, followed by an injection of resin to obtain the final composite part. To improve forming efficiency and avoid the "trial and error" process, numerical simulations are effective. The developed numerical approaches can be classified into kinematic and mechanical approaches [START_REF] Bussetta | Numerical forming of continuous fibre reinforced composite material: a review[END_REF][START_REF] Liang | A review of numerical analyses and experimental characterization methods for forming of textile reinforcements[END_REF][START_REF] Gereke | Experimental and computational composite textile reinforcement forming: a review[END_REF] The kinematic approaches are purely geometric and can give the fiber orientation during deformation [START_REF] Long | A simulation of reinforcement deformation during the production of preforms for liquid moulding processes[END_REF][START_REF] Wang | The draping of woven fabric preforms and prepregs for production of polymer composite components[END_REF] Although these methods are fast, they do not consider the mechanical behavior and therefore cannot simulate external process parameters such as draping force. Mechanical approaches [START_REF] Syerko | Models of mechanical properties/behavior of dry fibrous materials at various scales in bending and tension: a review[END_REF][START_REF] Baran | A Review on the Mechanical Modeling of Composite Manufacturing Processes[END_REF] have been developed at the different scales of textile composite materials. Compared with the method conducted at the microscale (on fibers) [START_REF] Durville | Simulation of the mechanical behaviour of woven fabrics at the scale of fibers[END_REF][START_REF] Latil | Towards the 3D in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres[END_REF] and mesoscale (on yarn) [START_REF] Lomov | Meso-FE modelling of textile composites: road map, data flow and algorithms[END_REF][START_REF] Nguyen | Mesoscopic scale analyses of textile composite reinforcement compaction[END_REF], macroscale modeling approaches consider the textile composites as a continuous medium and is thus more computationally effective. This continuity hypothesis has been verified in a variety of cases [START_REF] Gereke | Experimental and computational composite textile reinforcement forming: a review[END_REF][START_REF] Boisse | The difficulties in modeling the mechanical behavior of textile composite reinforcements with standard continuum mechanics of Cauchy. Some possible remedies[END_REF][START_REF] Khan | Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark[END_REF] The textile fabric can be considered as thin or thick, but generally, the thickness is much smaller than the two other dimensions. For this reason, shell finite elements can be adopted in macroscale approaches. Within the framework of classical shell elements (Kirchhoff or Mindlin), the bending rigidity is determined by the tensile stiffness and thickness. This cannot be applied to fibrous materials due to significant sliding occurring between the fibers, leading to a low bending rigidity. Such a weakness in bending has led to the development of membrane approaches (neglecting the bending behavior) [START_REF] Cao | An approach in modeling the temperature effect in thermo-stamping of woven composites[END_REF][START_REF] Lin | Predictive modelling for optimization of textile composite forming[END_REF][START_REF] Chen | Formability optimisation of fabric preforms by controlling material draw-in through in-plane constraints[END_REF] However, it has been shown that the bending behavior plays an important role in textile fabric forming, especially the formation and development of wrinkles [START_REF] Boisse | Simulation of wrinkling during textile composite reinforcement forming. Influence of tensile, in-plane shear and bending stiffnesses[END_REF][START_REF] Dangora | Predictive model for the detection of outof-plane defects formed during textile-composite manufacture[END_REF][START_REF] Thompson | Modelling defect formation in textiles during the double diaphragm forming process[END_REF][START_REF] Bai | Physics-based evaluation of the drapability of textile composite reinforcements[END_REF] To model this specific bending behavior, in the first proposed approach, a single-layer fabric was modeled as a superimposition of several elements [START_REF] Döbrich | Decoupling the bending behavior and the membrane properties of finite shell elements for a correct description of the mechanical behavior of textiles with a laminate formulation[END_REF][START_REF] Yu | A macroscale finite element approach for simulating the bending behaviour of biaxial fabrics[END_REF][START_REF] Dörr | Nonlinear hyperviscoelastic modelling of intra-ply deformation behaviour in finite element forming simulation of continuously fibre-reinforced thermoplastics[END_REF][START_REF] Dörr | A viscoelastic approach for modeling bending behavior in finite element forming simulation of continuously fiber reinforced composites[END_REF][START_REF] Haanappel | Formability analyses of uni-directional and textile reinforced thermoplastics[END_REF][START_REF] Jauffrès | Discrete mesoscopic modeling for the simulation of woven-fabric reinforcement forming[END_REF] The stress resultant shell approach [START_REF] Hamila | Semi-discrete shell finite elements for textile composite forming simulation[END_REF][START_REF] Chen | A hypoelastic stress resultant shell approach for simulations of textile composite reinforcement forming[END_REF] achieved the decoupling of bending and membrane behavior in a shell element. These approaches have been shown to be effective in predicting fabric mid-surface deformation, including the in-plane shear and out-of-plane bending deformation. Nevertheless, they are based on the classical shell theory which fails to predict the correct transverse shear strain for textile composites, as discussed later.

The transverse shear strain is the rotational displacement of material normals which are straight lines perpendicular to the initial mid-surface. The transverse shear strain is important as it will determine the strain https://doi.org/10.1016/j.ijmecsci.2023.108129 and stress state in the fibers at different thicknesses other than the midsurface. To predict the transverse shear strain, a direct approach would be to model multiple fiber layers as a stack of shell elements [START_REF] Soulat | Simulation of continuous fibre reinforced thermoplastic forming using a shell finite element with transverse stress[END_REF][START_REF] Huang | Simulation of wrinkling during bending of composite reinforcement laminates[END_REF] This approach becomes extremely time-consuming when the number of layers significantly increase. Different finite elements have been developed for this purpose, including solid elements [START_REF] Mathieu | Enhanced modeling of 3D composite preform deformations taking into account local fiber bending stiffness[END_REF][START_REF] Pazmino | Numerical modelling of forming of a non-crimp 3D orthogonal weave E-glass composite reinforcement[END_REF][START_REF] El Said | An iterative multiscale modelling approach for nonlinear analysis of 3D composites[END_REF] and fibrous shell elements beyond the classical shell kinematics [START_REF] Liang | A shell formulation for fibrous reinforcement forming simulations[END_REF][START_REF] Bai | A specific 3D shell approach for textile composite reinforcements under large deformation[END_REF] To study the structural response of laminated composites, such as buckling and vibration [START_REF] Kahya | Free vibrations of laminated composite beams with multiple edge cracks: numerical model and experimental validation[END_REF][START_REF] Youzera | Nonlinear damping and forced vibration analysis of laminated composite plates with composite viscoelastic core layer[END_REF][START_REF] Yan | Free vibration analysis of variable stiffness composite laminated beams and plates by novel hierarchical differential quadrature finite elements[END_REF], and thermomechanical properties [START_REF] Vidal | Thermo-mechanical analysis of laminated composite and sandwich beams based on a variables separation[END_REF][START_REF] Belbachir | Thermal flexural analysis of anti-symmetric cross-ply laminated plates using a four variable refined theory[END_REF][START_REF] Abualnour | Thermomechanical analysis of antisymmetric laminated reinforced composite plates using a new four variable trigonometric refined plate theory[END_REF], a lot of effort in recent years has been devoted to the development of solid-shell [START_REF] Xiong | A prismatic solid-shell finite element based on a DKT approach with efficient calculation of through the thickness deformation[END_REF][START_REF] Schwarze | A reduced integration solid-shell finite element based on the EAS and the ANS concept-Large deformation problems[END_REF][START_REF] Draiche | A computational investigation on flexural response of laminated composite plates using a simple quasi-3D HSDT[END_REF], single-layer (ESL) theories, zig-zag (ZZ) theories, Layerwise (LW) theories [START_REF] Liew | An overview of layerwise theories for composite laminates and structures: development, numerical implementation and application[END_REF][START_REF] Maji | Development and applications of shear deformation theories for laminated composite plates: an overview[END_REF][START_REF] Bakoura | A mechanical behavior of composite plates using a simple three variable refined plate theory[END_REF][START_REF] Allam | A generalized 4 -unknown refined theory for bending and free vibration analysis of laminated composite and sandwich plates and shells[END_REF], and Carrera unified formulation (CUF) models [START_REF] Yan | Free vibration analysis of variable stiffness composite laminated beams and plates by novel hierarchical differential quadrature finite elements[END_REF][START_REF] Wu | Geometrically nonlinear refined shell theories by Carrera Unified Formulation[END_REF] Most of these theories put effort into verifying the continuity of displacement (or stress) fields in the thickness. However, the specific behavior of dry reinforcements leads to significant sliding between plies. In order to take these effects into account, it is possible to use interface layers, but with a great increase in computational cost. High-order theories (second gradient solid elements [START_REF] Madeo | Thick fibrous composite reinforcements behave as special second-gradient materials: three-point bending of 3D interlocks[END_REF], Cosserat mediums [START_REF] Sakhaei | A finite deformation Cosserat continuum model for uncured carbon fibre composites[END_REF][START_REF] Dodwell | Internal wrinkling instabilities in layered media[END_REF]) have also been applied recently to dry reinforcements. However, it is still convenient and effective to apply classical shell theory, especially in commercial finite element codes.

The presented research had two objectives. The first was to present the difficulty of classical shell elements in predicting the transverse shear strain for fibrous composites. The four-point bending test of multilayer fibrous material was analyzed as direct evidence of these limitations. This is presented in Section 2. The second was to present a remedy approach within the framework of the classical shell element. In the proposed approach, the mid-surface was first determined by the stress resultant shell method (introduced in Section 3), then a complementary post-processing method was proposed in Section 4 to obtain the transverse shear. Section 5 compares numerical and experimental results in various bending and forming tests to show the effectiveness of the proposed approach.

The specificities of transverse shear behavior of fibrous reinforcement

In this section, the four-point bending test of fibrous reinforcement is presented and analyzed as direct evidence that the classical shell theory is not applicable for the transverse shear deformation of fibrous reinforcement.

The transverse shear in classical shell theory

Referring to an x 1 , x 2 , x 3 Cartesian coordinate system, x 1 and x 2 are on the shell mid-surface, and x 3 is the coordinate in the normal direction. In the classical shell elements, the normals are supposed to remain straight lines during the deformation. The bending curvatures χ and transverse shear strain Г are related by shell kinematics [START_REF] Auricchio | A triangular thick plate finite element with an exact thin limit[END_REF]:

χ = ⎡ ⎢ ⎢ ⎣ χ 11 χ 11 χ 12 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂θ 2 ∂x 1 - ∂θ 1 ∂x 2 ∂θ 2 ∂x 2 - ∂θ 1 ∂x 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Γ = [ γ 13 γ 23 ] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ θ 2 + ∂u 3 ∂x 1 -θ 1 + ∂u 3 ∂x 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (1)
Here θ 1 and θ 2 are the rotation of normals along the x 1 and x 2 directions, and u 3 is the mid-surface displacement in x 3 direction. Note that the curvature χ is not the mid-surface curvature χ which is the second derivative of the mid-surface displacement (χ 11 = -

∂ 2 u3 ∂ 2 x1 , χ 22 = ∂ 2 u3 ∂ 2 x2 , χ 12 = -2 ∂ 2 u3
∂x1∂x2 ). They are equal only under the condition that Γ = 0 (Kirchhoff shell). In the case of fibrous reinforcements, the bending stiffness is experimentally characterized by the mid-surface curvature and bending moment. Thus, Kirchhoff kinematics is implicitly used to decouple the bending behavior from the membrane behavior.

By adding transverse shear strain (Γ ∕ = 0), the Mindlin shell allows the shell normals to rotate and not always be perpendicular to the midsurface. However, the transverse shear strain γ α3 is determined by the transverse shear stiffness C s and the transverse shear stress resultant Q α3 (α takes value 1 or 2) [START_REF] Auricchio | A triangular thick plate finite element with an exact thin limit[END_REF][START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF][START_REF] Mindlin | Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates[END_REF]:

[ Q 13 Q 23 ] = C s [ γ 13 γ 23 ] (2)
It is shown next that the transverse shear stiffness C s does not apply to fibrous composite reinforcements. This was analyzed by four-point bending test.

The transverse shear of fibrous reinforcement in a four-point bending test

Fig. 1 shows the four-point bending test of fibrous materials (30 layers of paper in Figs. 1a, and 11 layers of interlock fabric G1151 in Fig. 1b). The two materials shared the same physical properties: the fibrous composition was almost inextensible, and relative slippage could occur between the adjacent layers. Both the mid-surface and material normals were marked to track the material deformation. It can be seen in Fig. 1 that the angle between the normal directors and mid-surface was far from 90 • Thus, it was straightforward to conclude that the Kirchhoff shell theory is not relevant. Moreover, in the Kirchhoff shell, only the mid-surface could fulfill the quasi-inextensibility conditions, and the material would be stretched/compressed on the top/bottom surface. This is also inconsistent with the physics of a fibrous material.

In a four-point bending test, the transverse shear stress resultant can be divided into three zones (neglecting gravity and friction). As shown in Fig. 2, in zone I, the sample was free without constraint, and the shear stress resultant was equal to zero. In zone II, the bending moment changed evenly, and there was a constant shear stress resultant. In zone III, the bending moment remained constant and the shear stress resultant stayed at zero (pure bending test).

Since the transverse shear stresses in zone I and zone III were zero, the shear strain given by the Mindlin shell would according to Eq. ( 2) also be zero. This differed from the experimental results (Fig. 1). A comparison of the shear strain is given in Table 1. The inconsistent shear strain by the Mindlin shell theory did not ensure the quasiinextensibility of the fibers at different thicknesses. From the fourpoint bending test, it was found that the transverse shear strain in the fibrous materials should not be governed by the constitutive equation (Eq. ( 2)), but rather by the quasi-inextensibility of the fibers and the slippage between them.

Numerical result of the four-point bending test with the Mindlin shell element

To corroborate the theoretical analysis, numerical simulations of four-point bending tests were carried out in the commercial software Abaqus. In the simulation model, the material thickness was 15 mm, and the length was 300 mm, i.e., the same as in the experiment of Fig. 1b. A Mindlin shell element (S4) was used, and the element size was set as 5 mm × 5 mm. In order to take into account the friction between fibers, the bending stiffness of 11 layers of G1151 was adopted in the simulation measured by a three-point bending test [START_REF] Bai | A specific 3D shell approach for textile composite reinforcements under large deformation[END_REF] The transverse shear stiffness in the (x 1 , x 3 ) plane was set from 0.001 N/mm to 100 N/mm. Since the shell transverse shear strain could not be directly visualized in Abaqus, the shell cross-section was plotted in Matlab using the numerical result of the shell mid-surface position and the orientation of material normals, both given by Abaqus. The simulation is presented in Fig. 3. Regardless of the value of the transverse shear stiffness, the transverse shear strains in Zone I and III were zero, which corresponded to the theoretical results in the Mindlin shells.

As depicted in Fig. 3a, the transverse shear stiffness was almost equal to zero, and the normals were found to be quasi-vertical. In this situation, the orientation of the normals was close to the fibrous material shown in Fig. 1. However, the mid-surface deformation was far from the

Table 1

Comparison of the transverse shear strain γ13 in the four-point bending test. experimental shape, and the reason for this was that the bending deformation was replaced by transverse shear due to the low transverse shear stiffness. Moreover, the transverse shear strain determined the inplane tensile strain at different thicknesses. According to Fig. 4, the quasi-inextensibility of the fibers would not be ensured due to the incorrect transverse shear strain. The simulations confirm that neither Kirchhoff nor Mindlin is suitable for describing the transverse shear of a stack of fibrous material. Nevertheless, the classical shell theory was quite effective for modeling the specific in-plane shear and bending behavior to obtain the midsurface deformation. The following section proposes a post-processing method to correct the transverse shear given by stress resultant shell simulations that lie within the Kirchhoff theory.

The stress resultant shell approach for the mid-surface deformation

The stress resultant shell approach was used to simulate the midsurface deformation by considering the independent bending behavior of the textile reinforcement. With this aim, the internal virtual work δW int was regarded as the sum of virtual works of the in-plane ten-sionδW tension , in-plane shear δW shear and out-of-plane bendingδW bending :

Here, A is the shell mid-surface. δε 11 , δε 22 are the axial strain increments, δε 12 is the in-plane shear strain increment, and δχ 11 , δχ 22 ,δχ 12 are the bending curvature increments. N 11 , N 12 , N 22 are the stress resultants and M 11 , M 12 , M 22 are the stress moments for which the direction is shown in Fig. 5. The coupling component M 12 was neglected due to the bending stiffness mainly originating from the warp and weft yarns. Hypoelastic laws were adopted in the simulation approach to compute the objective time derivative of the stress resultants N ∇ and the stress momentsM ∇ :

N ∇ = C : ε M ∇ = D : χ (4)
With Eq. ( 4), the stress resultant N was linked to the membrane strain ε through the membrane stiffness matrix C, and the stress moment M was linked to the bending curvature χ through the bending stiffness matrix D. The membrane and bending response were thus naturally independent thereby satisfying the requirement of the textile material.

The stress resultant shell approach has been implemented in the Abaqus/Explicit subroutine VUGENS [START_REF] Chen | A hypoelastic stress resultant shell approach for simulations of textile composite reinforcement forming[END_REF] According to the internal virtual work expressed in Eq. ( 3), no transverse shear energy is considered in the simulation approach. However, to account for the friction between fibers during transverse shearing, the bending stiffness used in Eq. ( 4) is measured on a multilayer material, rather than bending deduced from a single-layer experiment. This is a simple and effective way to obtain the bending deformation of the mid-surface without complicating the approach [START_REF] Liang | A shell formulation for fibrous reinforcement forming simulations[END_REF].

To verify the normal transformation in the stress resultant shell approach, one can conduct a numerical simulation of the four-point bending test. The simulation settings were the same as those given in Section 2.3, and the simulation result can be seen in Fig. 6. The midsurface bending deformation was close to the experiment (shown in Section 5.2.2). However, it was found that the normals remained perpendicular to the deformed mid-surface which was not consistent with the fibrous material shown in Fig. 1.

The kinematic modeling approach for transverse shear of textile fabric

After determining the mid-surface position by the stress resultant shell method, the transverse shear would be geometrically calculated based on the quasi-inextensibility of the fibers. This is presented in the next section.

Kinematic hypothesis

The fibrous domain that contains parallel fibers is considered in this research. The objective of the modeling approach was to re-calculate the transverse shear strain in the fibrous reinforcement draping. Three kinematic hypotheses were put forward according to results from bending experiments (Fig. 1):

• The material director remained in straight lines during the fabric deformation, but was not necessarily perpendicular to the midsurface. The hypothesis of straightness was not fulfilled when the reinforcement became very thick, but it could be verified for low or medium thicknesses which are representative of most cases involving fibrous reinforcement drapings.

• The thickness along the material director could be stretched. However, neglecting the delamination in some draping cases, the material thickness h in the normal direction remained constant as it was the sum of the thicknesses of thousands of fibers. The obtained material normals of the nodes in the mid-surface are marked blue.

• Parallel fibers could slip in relation to each other along the thickness, but the in-plane slippage between the warp and weft yarn directions was neglected.

Modeling process

The numerical simulation of the mid-surface gave the position of the nodes in the mid-surface and the orientation of the normals under the global coordinate system. According to the previous analysis, the output normals remained perpendicular to the deformed mid-surface. Thus, the position of the nodes in the bottom and top surfaces can first be determined with:

V ( B k m ) = V(B) + ( -1) α n B h 2‖ n B ‖ (5)
Here, V(B k m ) is the position vector of the node B k m , m belongs to the set (1, 2) (m = 1 represents the bottom surface and m = 2 represents the top surface), V(B) is the position vector of a random node B in the midsurface, h is the shell thickness which remains constant with the initial value and n B is the shell normal vector output from the simulation. As shown in Fig. 7a, the shell normals n A ,n B are perpendicular to the midsurface, and the nodes

A k 1 , A k 2 , B k 1 , B k
2 calculated by the Eq. ( 5) were thus marked as Kirchhoff nodes with the superscript k in this research. To ensure the quasi-inextensibility of the fibers in the top and bottom surfaces, they had to be adjusted to the fibrous nodes (marked with the superscript f).

In the case of drapings in one fiber direction, the Kirchhoff node B k m , was adjusted along the fiber direction to the fibrous node B f m according to (shown in Fig. 7a):

( x A k m -x B f m ) 2 + ( y A k m -y B f m ) 2 = L 2 AB K ( x B f α , y B f α ) = 0 (6)
Here, A, B are the nodes in the mid-surface, and A k m is the Kirchhoff node whose normals remained perpendicular to the deformed midsurface under the applied boundary conditions. The adjustment needs to start from these starting points to ensure that the calculated fibrous nodes in the top and bottom surfaces corresponded to each other.

K(x B f α , y B f α ) is the line equation of A k m B k
m determined by Eq. ( 5), and L AB is the distance between points A and B. Since the mid-surface deformation was simulated by the classical shell elements with a small length, the straight line distance was used in Eq. ( 6) for the sake of simplicity.

In the case of drapings in two fiber directions, the adjustment of the Kirchhoff nodes occurs simultaneously along the weft and warp yarn directions. In order to obtain the starting points for each yarn, the adjustment is first conducted on the points that lay on a given boundary condition. For these points, they only need to be adjusted along one fiber direction. As shown in Fig. 7b, they are points B and C. Referring to point A, whose normal keeps perpendicular to the deformed mid-surface, Eq.

( x B f m -x D f m ) 2 + ( y B f m -y D f m ) 2 + ( z B f m -z D f m ) 2 = L 2 BD ( x C f m -x D f m ) 2 + ( y C f m -y D f m ) 2 + ( z C f m -z D f m ) 2 = L 2 CD F ( x D f m , y D f m , z D f m ) = 0 (7)
Here, L BD and L CD are the line distance.

F(x D f m , y D f m , z D f m
) is the top or bottom surface function determined by the nodes obtained in Eq. ( 5). Based on the inextensibility of the fibers, Eq. ( 7) has been adopted in the kinematic approaches for the calculation of the mid-surface deformation [START_REF] Long | A simulation of reinforcement deformation during the production of preforms for liquid moulding processes[END_REF][START_REF] Wang | The draping of woven fabric preforms and prepregs for production of polymer composite components[END_REF] After determining the fibrous nodes D f m , the orientation of the normals was obtained by connecting the corresponding fibrous nodes in the bottom and top surfaces.

The computation process is purely geometric, giving this approach a low computation cost. This approach combines the advantage of mechanical approaches and kinematic approaches. It can also be an extension of other methods [START_REF] Döbrich | Decoupling the bending behavior and the membrane properties of finite shell elements for a correct description of the mechanical behavior of textiles with a laminate formulation[END_REF][START_REF] Yu | A macroscale finite element approach for simulating the bending behaviour of biaxial fabrics[END_REF][START_REF] Dörr | Nonlinear hyperviscoelastic modelling of intra-ply deformation behaviour in finite element forming simulation of continuously fibre-reinforced thermoplastics[END_REF][START_REF] Dörr | A viscoelastic approach for modeling bending behavior in finite element forming simulation of continuously fiber reinforced composites[END_REF][START_REF] Haanappel | Formability analyses of uni-directional and textile reinforced thermoplastics[END_REF][START_REF] Jauffrès | Discrete mesoscopic modeling for the simulation of woven-fabric reinforcement forming[END_REF] developed to obtain mid-surface deformation. Based on the kinematics mechanism, the approach extends the deformation in the mid-surface to the planes in the different thickness positions, which can be meaningful for the simulation of fabric forming with significant nonzero transverse shear strain.

Comparison of numerical and experimental results

The objective of the proposed modeling approach was to simulate the transverse shear deformation in the textile composite draping. In order to verify its correctness and effectiveness, several bending and hemispherical forming tests were carried out and compared with the numerical results obtained by the kinematic modeling approach.

Introduction to the studied textile

The textile reinforcement studied in this work was the interlock carbon fabric G1151. The geometry parameters are listed in Table 2.

Table 2

Geometry parameters of interlock fabric G1151. The properties of a single-ply G1151 fabric are listed in Table 3. The in-plane shear behavior was obtained through the bias-extension test [START_REF] Lebrun | Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics[END_REF][START_REF] Harrison | Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Cao | Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results[END_REF][START_REF] Schirmaier | Characterisation of the draping behaviour of unidirectional non-crimp fabrics (UD-NCF)[END_REF], and the bending behavior was measured from the cantilever test under gravity [START_REF] Syerko | Models of mechanical properties/behavior of dry fibrous materials at various scales in bending and tension: a review[END_REF][START_REF] Chen | Analysis of the forming of interlock textile composites using a hypoelastic approach[END_REF][START_REF] De Bilbao | Experimental study of bending behaviour of reinforcements[END_REF][START_REF] Sachs | Viscoelastic bending model for continuous fiber-reinforced thermoplastic composites in melt[END_REF] A very large and constant tensile stiffness was chosen as input to ensure the quasi-inextensibility of the fibers during the numerical simulation.

Bending tests of fibrous stacks

In this section, the specimens consist of parallel fibers (paper or reinforcement) that can slide relative to each other. Cantilever and fourpoint bending tests with large displacement are considered in order to obtain significant transverse shear deformation. Experimental results are then compared to numerical simulations.

Cantilever bending test

The cantilever bending test was conducted on a book consisting of 130 pages [START_REF] Liang | A shell formulation for fibrous reinforcement forming simulations[END_REF] The properties of a single sheet of paper are given in Table 4. Since the friction between the sheets was negligible, the properties of the stack could be deduced from it. The experimental geometry is given in Fig. 8a. The specimen is bent under gravity and the experimental bending result is shown in Fig. 8b. The mid-surface simulation was conducted using the stress resultant shell approach. The simulation setting was the same as that described in Section 2.3. The simulation result of the mid-surface is presented in Fig. 9a, and was found to be close to the experimental Subsequently, the directors of the normals were computed using the modeling approach presented in Section 4, and the obtained cross- section is shown in Fig. 8c. Four positions (Fig. 8b) were selected to compare the rotation angle and thickness of normals between experimental and numerical results. The compared results of the rotation angle, given in Fig. 9b, showed good agreement. The thickness of the normals, cf Fig. 9c, exhibited the same tendency with some deviations.

The experimental thickness of the normals was higher than those obtained the by modeling approach. This may be related to delamination occurring during the experiment.

Four-point bending test

The four-point bending test was performed on a stack of 11 layers of interlock fabric G1151 in which the fibers in each layer were parallel to each other. The experimental geometry is given in Fig. 10a. The punch displacement was 50 mm, and the final deformed shape is shown in Fig. 10b.

Fig. 11a displays a comparison between experimental and numerical results of the mid-surface shape. Good agreement was found. The calculated cross-section through the modeling approach is shown in Fig. 10c. The rotation angle and thickness of the normals were compared with experimental results in Fig. 11b and11c, respectively, and were found to be close, which cannot be obtained by classical shell elements.

Hemisphere forming experiment

The hemisphere is a benchmark shape for the textile reinforcements forming analysis [START_REF] Ten Thije | A multi-layer triangular membrane finite element for the forming simulation of laminated composites[END_REF][START_REF] Peng | Validation of a non-orthogonal constitutive model for woven composite fabrics via hemispherical stamping simulation[END_REF][START_REF] Khan | A parametric sensitivity study on preforming simulations of woven composites using a hypoelastic computational model[END_REF][START_REF] Liu | Experimental characterization and analysis of fiber orientations in hemispherical thermostamping for unidirectional thermoplastic composites[END_REF][START_REF] Aridhi | Textile composite structural analysis taking into account the forming process[END_REF][START_REF] Bai | The fibrous shell approach for the simulation of composite draping with a relevant orientation of the normals[END_REF][START_REF] Chen | A hypoelastic approach for simulating the independent bending behavior of textile composite within the stress resultant shell[END_REF]. The experimental setting is given in Fig. 12. The specimen was placed between the blank holder and the die, the final displacement of the punch was set at 75 mm, and the forming speed was 2.5 mm/s. To observe and record the material deformation, the blank holder and die were made transparent and two cameras were positioned: one at the top and the other at the bottom.

Two types of fabric stack sequences were investigated separately in the experiment, including a four-layer stack in 0 • /90 • configuration and a four-layer stack in ±45 • configuration. The final deformed fabric shapes are portrayed on the right-hand sides of Figs. 14a and15a. None of them exhibited wrinkles.

12. Schematic diagram of the forming experiment. The material is placed between transparent blank holder and die mold. Two cameras are used to capture the material deformation on the top and bottom surfaces. The corresponding simulation to obtain the mid-surface was carried out as follows. The blank holder, die, and punch were modeled as rigid bodies. Triangle shell elements were adopted to mesh the fabric, and the element side length was 5 mm × 5 mm. The Coulomb friction model was used to describe the contact between the fabric and the forming tools, and the coefficient was 0.2. The deformation on the top and bottom surfaces was calculated using the proposed modeling approach, and the numerical shape of the four-layer ±45 • fabric is presented at the bottom of Fig. 13.

Comparison of in-plane shear angle

The experimental shear angle was calculated according to the painted marker points, and the numerical shear angle at the top and bottom surfaces were computed from adjacent fibrous nodes. From the experimental results (top of Fig. 13), the in-plane shear angle at the top and bottom surfaces of the fabric differed somewhat even when the material remained flat between the blank holder and die mold. This different shear angle would also cause a rotation of the material normals. The presented modeling approach was able to reproduce the inplane shear variation along the thickness. Numerical tendencies were the same as their experimental counterparts, as shown in Fig. 13.

Comparison of rotation of normals

To study the rotation of the normals in the hemispherical forming, eight marker points located in the fiber direction and the symmetry direction were chosen. The rotation angle of the normals in the experiment was calculated using the corresponding marker points on the top and bottom surfaces. The two different stack sequences were then analyzed.

(a) Four-layer stack at 0 • /90 • The selected marker points (labeled 1 to 8) in the hemispherical experiment of the four-layer stack at 0 • /90 • are shown in Fig. 14a. The hemispherical punch had a centrosymmetric shape, and the initial shape of the specimen was a 300 mm × 300 mm square. The position of the marker points thus changed during deformation, whereas the crosssection of the specimen remained in the vertical plane.

The numerical results of the cross-section along the fiber and symmetry directions are presented in Fig. 14b and14d. The normal directors at the selected marker points (No.1-8) are also indicated. Point No.1 remained between the die and blank holder after the deformation. The material normals stayed vertical, which was a consequence of the specific punch shape and the inextensibility of the fibers. Points No.2-4 were located in the hemispherical zone after deformation. The rotation of the normals was found to be significant which was essential to ensure the quasi-inextensible of the fibers. At points No.5-8, there was no inextensibility constraint because of the absence of fibers. The rotation of the normals occurred due to the shear deformation being different between warp and weft yarns at the top and bottom surfaces.

As shown in Fig. 14c and14e, the numerical results of the rotation angle were found close to the corresponding experimental results. The direction of material normals was important as it determined the strain of the fibers in the different thickness positions. Unlike the proposed kinematic modeling approach that can correctly predict the direction of normals in hemisphere forming, other approaches listed in the introduction only focus on the mid-surface deformation. Based on the classical shell kinematics, these alternative approaches would lead to the normals being perpendicular to the mid-surface after the deformation which would be incorrect for a fibrous medium.

(b) Four-layer stack at ±45 • In the hemispherical experiment of the four-layer stack at ±45 • , Fig. 15a shows the location of the selected marker point. Points No.1-4 were in the symmetry direction while points No.5-8 were placed along the fiber direction. The numerical results of the cross-section where the marked points were located are illustrated in Fig. 15b and15d. The orientation of the normals obtained by the modeling method is presented in Fig. 15c and15e, and as can be seen, the outcome was in fair agreement with the experimental results.

Discussion

A method was presented to deduce the response of multi-layer fabrics from the mid-surface solution. After calculating the mid-surface deformation through finite element analysis, it was in certain cases possible to compute transverse shear using a kinematic approach relying on the inextensibility of the fibers. This approach decoupled the calculation of transverse shear from the bending curvature of the mid-surface. For simple forming shapes, the solution of this two-step calculation was the same as the more complex approaches like 3D shell elements [START_REF] Bai | The fibrous shell approach for the simulation of composite draping with a relevant orientation of the normals[END_REF]. Nevertheless, it should be noted the limitation of the presented approach.

As introduced in Section 4.2, the presented algorithm requires the normals rotation of the starting point and the two fiber directions crossing the starting point as boundary conditions (BCs). For a case of shapes with two symmetry planes (as in hemispherical forming), these BCs are given: the normal rotation of the central node is zero, and the two fiber directions are determined by the symmetry planes. In 3D forming cases with less than two symmetry planes, the BCs are not explicitly given, and additional assumptions have to be made. Note that this is not a specific limitation of the presented approach, but it is common to other kinematic approaches for calculating the mid-surface deformation [START_REF] Pickett | Simplified and advanced simulation methods for prediction of fabric draping[END_REF][START_REF] Hancock | The use of kinematic drape modelling to inform the hand lay-up of complex composite components using woven reinforcements[END_REF].

Regarding normals rotation of the starting point, which is specific to our multi-layer approach, further comments can be made. Let us suppose that this BC is not given by symmetry conditions, and has to be assumed. In this case, the algorithm will provide a family of solutions depending on the chosen BC. As illustrated in Appendix A, the prescribed BC will affect all the material directors in a non-linear way, but all solutions can be mathematically related to one another based on the quasi-inextensibility of the fibers. The final choice of the correct solution could depend on an energy calculation. For example, the assumption can be made that the best solution is the one that minimizes the integral of transverse shear over the fabric. This is somehow related to the minimization of dissipated energy from slippage between the plies of fabric and the forming tools.

Conclusions

Through the example of a four-point bending test, this work first presented the limitation of classical shell theories for forming simulations of textile composite reinforcements. The Kirchhoff or Mindlin shell theories can be applied to simulate the mid-surface deformation, but neither of these theories can provide the correct transverse shear in fibrous mediums. It is quite straightforward to establish that the Kirchhoff shell was not relevant because of the zero transverse shear hypothesis. Although the Mindlin theory authorizes transverse shear in shells, it is not able to model transverse shear in fibrous stacks due to the fact that the transverse shear should be governed by fiber inextensibility in each layer outside the mid-plane, rather than by transverse shear stiffness used in Mindlin shells. The transverse shear deformation is important as it will determine the strain of fibers in the different thickness positions.

In order to take advantage of the classical shell elements to describe the specific bending behavior of fibrous material, a kinematic modeling approach was proposed to calculate the transverse shear deformation without complicating the shell approach. After simulating the midsurface response by stress-resultant finite element shell, the top and bottom layer (and layer-wise) response of the fabric was deduced from fiber inextensibility. This approach made it possible to extend the deformation in the mid-surface to planes of different thicknesses. The computation process is purely geometric, giving this approach a low computation cost. Experimental bending tests and hemispherical forming were compared to the solutions of the proposed approach, with good agreement, proving its correctness and effectiveness. The modeling approach can be implemented as a post-processing application, and can therefore be used by any user, together with any finite element software.

Appendix A. The solutions under varying normal rotations of the starting point

A 2D example, with a unique already known fiber direction, is illustrated here. With Eq. (A.4), it is shown that all solutions can be mathematically related to one another based on the quasi-inextensibility of the fibers.

Fig. 1 .

 1 Fig. 1. Transverse shear in the four-point bending test of fibrous material. Two materials that shared the same bending physics are used: (a) 30 layers of paper. (b) 11 layers of interlock fabric G1151. In order to observe the no-zero transverse shear deformation, the mid-surface and material normals are all marked in the test.

Fig. 2 .

 2 Fig. 2. Transverse shear stress resultant distribution in the four-point bending test. The shear stress is divided into three zones I, II, III. In zone I and II, the shear stress is zero, where constant in zone III.

Fig. 3 .

 3 Fig. 3. The four-point bending simulation using Mindlin shell element with different transverse shell stiffness. (a) Cs = 0.001 N/mm (b) Cs = 0.005 N/mm (c) Cs = 0.025 N/mm (d) Cs = 100 N/mm. In these four tests, all the simulation settings and other material properties are the same. To visualize the transverse shear, the simulation results are post-processed in Matlab to plot the material normals.

Fig. 4 .

 4 Fig. 4. Variation of the in-plane tensile strain over thickness in the Mindlin shell. The simulation result corresponds to Fig. 3d, where a relatively large transverse shear stiffness is used. The color bar represents tensile logarithmic strain.

Fig. 5 .

 5 Fig. 5. The elementary cell of woven reinforcement. Arrows represent the direction of stress resultants N 11 , N 22 , N 12 and moments M 11 , M 22 in the woven reinforcement that is interlaced by the warp and weft yarn. These components rotate along with the yarn direction under in-plane shear deformation.

Fig. 6 .

 6 Fig.6. The orientation of the normals in the four-point bending test using the stress resultant shell approach. The angle between material normals and mid-surface is highlighted by red arrows at two different positions.

Fig. 7 .

 7 Fig. 7. Schematic diagram of the calculation of the fibrous nodes in the top and bottom surfaces (a) Draping in one fiber direction. (b) Draping in two fiber directions. The obtained material normals of the nodes in the mid-surface are marked blue.

( 6 )

 6 is used to obtain the corresponding fibrous nodes B f m and C f m . For a random point D, the adjustment to obtain the coordinate of the fibrous nodes D f m is done with the following equation:

  C αα = 21500 N/mm dN αα = C αα dε αα (α= 1 or 2) In-plane shear behavior dN11 = C12dγ C12 = K0 + K1|γ| + K2|γ| 2 + K3|γ| 3 + K4|γ| 4 K0 = 0.09 N/mm K1 = -0.83 N/mm K2 = 2.92 N/mm K3 = -4.01 N/mm K4 = 2.03 N/mm Bending behavior D αα = 7.45 Nmm dM αα = D αα dχ αα (α= 1 or 2)

Fig. 8 .Fig. 9 .

 89 Fig. 8. 130-layer paper cantilever bending test under gravity. (a) Geometry. (b) Experimental results[START_REF] Liang | A shell formulation for fibrous reinforcement forming simulations[END_REF]. In the experiment, the material normals that are perpendicular to the initial mid-surface are marked black. (c) Numerical results using resultant shell finite element calculation, and kinematic approach for the calculation of transverse shear. Four different positions 1-4 are marked in order to compare the material normals.

Fig. 10 .

 10 Fig. 10. Four-point bending of an 11-layer interlock fabric G1151. (a) Geometry. (b) Experimental results. In the experiment, the mid-surface and material normals are marked white. (c) Numerical results using resultant shell finite element calculation, and kinematic for the calculation of transverse shear. Five different positions 1-5 are marked to verify the correctness of numerical results.

Fig. 11 .

 11 Fig. 11. Comparisons between experimental and numerical results in four-point bending test of an 11-layer interlock fabric G1151. (a) The mid-line. (b) The rotation of the normals. (c) The thickness.

Fig. 13 .

 13 Fig. 13. Experimental and numerical results of the shear angle in hemispherical forming of a four-layer ±45 • fabric. The upper part of the figure shows the experimental results: left-hand image is given by the bottom camera (bottom view), right-hand image is given by the top camera (top view). The lower part of the figure gives the corresponding numerical results on top and bottom layers. Shear angle γ is compared at two different positions.

Fig. 14 .

 14 Fig. 14. Four-layer 0 • /90 • hemispherical forming. (a) Initial and deformed fabric, the selected eight marker points are located in the fiber and symmetry directions. (b)(d) Numerical results of the cross-section along the fiber and symmetry directions. The orientation of normals at the selected maker points is highlighted with red lines. (c)(e) Normals rotation comparison along the fiber and symmetry directions.

Fig. 15 .

 15 Fig. 15. Four-layer ±45 • hemispherical forming. (a) Initial and deformed fabric, the selected eight marker points are located in the symmetry and fiber directions. (b)(d) Numerical results of the cross-section along the symmetry direction and fiber direction. The orientation of normals at the selected maker points is highlighted with red lines. (c)(e) Normal rotations comparison along the symmetry direction and fiber direction.

Fig. A. 1 .

 1 Fig. A.1. Simulation of cantilever bending test with three different prescribed rotations at the clamped edge. (a) 0 • rotation. (b). 15 • rotation. (c) 30 • rotation. The dashed line in the figure represents the vertical direction. Four positions 1-4 are selected to compare to the rotation of normals under different boundary conditions.

Fig. A. 2 .

 2 Fig. A.2. Normal rotations under three different boundary conditions of cantilever bending test. The double arrow shows the difference in rotation angle at different node positions.

  Fig A.1 presents three different solutions of a cantilever test with the same mid-plane solution given by the finite element analysis but three different normal rotations β at the clamped edge. The normal rotations at other positions are simulated using the proposed modeling approach. It is shown in Fig A.2 that the prescribed BC will affect all the material directors in a non-linear way.After determining the corresponding node positions [x Top ,y Top ][x Bot ,y Bot ] on the top and bottom surfaces, respectively, the normal orientation θ can be calculated according to:θ = arctan ( y Topy Bot x Topx Bot ) (A.1)When there is a prescribed normal rotation β at the starting point, all node positions are updated based on fiber inextensibility. Displacement of the top/bottom nodes are marked as L c . Top , α Bot are the top and bottom angles between the fiber direction and the horizontal line, and can be obtained by the adjacent nodes on the top or bottom surfaces. By combining Eq A.1-3, the new orientation of normal θ c can be determined by:θ c = arctan ⎛ ⎜ ⎝ y Top -y Bot -h 2 cotβ ( sinα Top + sinα Bot ) x Topx Bot + h 2 cotβ ( cosα Top + cosα Bot )

Table 3

 3 Mechanical properties of the single layer fabric G1151.

Table 4

 4 Properties of a single sheet of paper.

	Properties Tensile stiffness Bending stiffness Areal density Thickness	Values C = 230 N/mm D = 0.4 N mm 81 g/m 2 0.1 mm

Acknowledgment

This work was supported by the ANR, (French Research Agency), grant ANR-18-CE06-0011-04-AMOC, and by the China Scholarship Council (CSC) (B. Chen).

Data availability

Data will be made available on request.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.