# A hyperelastic approach for modeling the membrane behavior in finite element forming simulation of unidirectional non-crimp fabrics (UD-NCF)

Bastian Schaefer, Dominik Doerr, Ruochen Zheng, Naïm Naouar, Luise Kaerger

## ► To cite this version:

Bastian Schaefer, Dominik Doerr, Ruochen Zheng, Naïm Naouar, Luise Kaerger. A hyperelastic approach for modeling the membrane behavior in finite element forming simulation of unidirectional non-crimp fabrics (UD-NCF). Composites Part A: Applied Science and Manufacturing, 2024, 185, 10.1016/j.compositesa.2024.108359. hal-04706377

HAL Id: hal-04706377

https://hal.science/hal-04706377v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A hyperelastic approach for modeling the membrane behavior in finite element forming simulation of unidirectional non-crimp fabrics (UD-NCF)

Bastian Schäfer<sup>a,*</sup>, Dominik Dörr<sup>b</sup>, Ruochen Zheng<sup>c</sup>, Naim Naouar<sup>c</sup>, Luise Kärger<sup>a</sup>

## ARTICLE INFO

### Keywords:

A. Fabrics/textiles

A. Unidirectional non-crimp fabric

C. Process simulation

E. Forming

E. Preforming

## ABSTRACT

Unidirectional non-crimp fabrics (UD-NCFs) are highly suitable for high-performance components due to their excellent lightweight potential. However, during forming they are prone to wrinkling and gapping compared to woven or biaxial textiles. Macroscopic simulation models can be used to efficiently predict these effects as well as the global forming behavior for complex geometries. Therefore, a new hyperelastic membrane model is proposed to describe the typical deformation of UD-NCFs based on superimposed shear, transverse tension and compression perpendicular to the fiber rovings. The model is parameterized using the forces and different ratios of superimposed strains obtained in experimental off-axis-tension-tests at bias angles of 30°, 45° and 60°. The resulting approach is validated by forming simulations of a hemisphere and tetrahedron geometry in different configurations and quantitative comparison to experimental tests. The model accurately predicts the forming behavior of UD-NCF with a good agreement of the global deformation behavior and local strains.

## 1. Introduction

Unidirectional non-crimp fabrics (UD-NCFs) have straight fibers without undulations common in woven fabrics and therefore provide a higher lightweight potential. They consist of one layer of unidirectional fiber rovings linked together with polymer stitching. The lack of a second reinforcement direction for stabilization impedes handleability, and producing good forming results is difficult. UD-NCFs are prone to wrinkling, gapping or fiber waviness during forming of complex geometries [1–5]. Thus, macroscopic modeling approaches are used to predict these defects and the global forming behavior efficiently [6–9].

Woven fabrics have often been the focus of research on macroscopic forming simulations due to their better formability, as evident from the comprehensive reviews in literature [10–15]. Thereby, hypoelastic or hyperelastic approaches are most commonly used to model their membrane behavior. Hypoelastic approaches require a non-orthogonal constitutive model to account for fiber rotation and typically relate the stress rate to the strain rate directly via a non-linear stiffness [16–18]. Hyperelastic approaches define a strain energy density as a potential for the stress, which can be additively decomposed based on directional strain-invariants representing each deformation mode [19–22]. Additionally, biaxial tension as well as tension-shear interactions were investigated for woven fabrics by introducing suitable couplings between the deformation modes [23–26].

In contrast, biaxial NCF [27–31] or UD-NCF [6–8,32–35] have been investigated much less and often based on mesoscopic approaches [32–34]. Mesoscopic approaches discretely model the material's constituents and can only be used to a limited extent for component forming simulations due to their high numerical effort. Nevertheless, they can predict local forming effects, while macroscopic approaches capture the deformations mechanisms occurring on the meso-scale in a homogenized manner to indicate areas with a high likelihood of e.g. waviness, gapping or in-plane compaction [7].

Membrane approaches for NCFs need to consider not only shear, which is the main deformation mode for woven fabrics during forming [13], but also stitching deformation and roving slippage. Yu et al. [36] developed a non-orthogonal hypoelastic model for Biax-NCF considering an asymmetrical shear behavior due to the tricot stitching oriented at 45° relative to the fibers. They highlight the influences of friction between the blanks and tools, blankholder forces and blank shape in numerical forming studies on a hemispherical shape. Similar approaches were later developed by Chen et al. [37] and Guzman-Maldonado et al. [38]. Khiêm et al. [30] proposed a hyperelastic model as a user-defined material in ABAQUS for a Biax-NCF with a chain stitch oriented at 45°. The model considers stitching and fiber tension as well as in-plane shear. Fiber slippage was neglected in the approach, but

<footnote><sup>a</sup> Lightweight Engineering, Institute of Vehicle System Technology (FAST), Karlsruhe Institute of Technology (KIT), Karlsruhe, 76131, Germany</footnote>

<footnote><sup>b</sup> Simutence GmbH, Karlsruhe, 76131, Germany</footnote>

<footnote><sup>c</sup> Université de Lyon, LaMCoS, CNRS, INSA Lyon, F-69621, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: bastian.schaefer@kit.edu (B. Schäfer), luise.kaerger@kit.edu (L. Kärger).</footnote>
despite the simplicity of the approach good agreement with experimental forming tests on a double dome geometry was achieved. Mallach et al. [31] modeled the behavior of a Biax-NCF primarily based on a nonlinear shear stress-strain curve in PAM-FORM. Acceptable results were obtained for the global forming behavior of a complex geometry for different blankholder configurations, but neglecting roving slippage resulted in noticeable deviations in the fiber orientation in areas of large shear deformation. Bel et al. [39] modeled inter-ply sliding between the two fiber directions of a Biax-NCF with a semi-discrete approach, utilizing two element layers connected with bar elements representing the stitching and enabling sliding in the fiber direction.

All approaches for Biax-NCFs neglect intra-ply slippage between rovings and large tensile strains of the stitching, which are highly relevant mechanisms for UD-NCF during forming due to the lack of a second fiber direction [1,40,41]. Therefore, the exclusive focus on the shear behavior during characterization and modeling is not sufficient. This is also evident in the study of Krogh et al. [35], in which the ABAQUS built-in FABRIC material model is utilized to describe the shear behavior of a quasi-unidirectional glass-fiber NCF. The model predicts the early force-displacement-relation well in a 45°-off-axis-tension test (OAT), also called bias-extension-test. However, the model overestimates the shear response for larger deformations as UD-NCF deforms more under simple shear instead of pure shear. A very comprehensive macroscopic approach was proposed by Schirmaier et al. [6]. It consists of an elastic-plastic transverse tensile behavior coupled to the elastic-plastic shear behavior via a 2D-yield surface and is superimposed with a 2D-elastic compressive stiffness perpendicular to the fiber rovings which is also coupled to the shear angle. The model is inversely parameterized using 30°-, 45°-, and 60°-OATs, approximating the applied forces well but underestimating the local strains in the main deformation zones of the specimen centers [8]. Nevertheless, it achieves good results in qualitative [6] and quantitative [7] comparisons to component forming results. However, the complex model requires a high number of material parameters to describe the membrane behavior, resulting in a difficult and time-consuming parameterization process and limited applicability to other materials. Schäfer et al. [8] proposed a simplified hyperelastic approach based on similar principles as the model of Schirmaier et al. [6], but omitted the couplings and prescribing only non-linear elastic stiffnesses. It was used to investigate the minimal level of complexity necessary to describe the most relevant deformation mechanisms of UD-NCF. Thereby, it was shown that at least a coupled transverse tensile behavior is necessary to approximate the large stitching deformation and relative fiber slippage in a homogenized way. Ghazimoradi and Montesano [9] recently parameterized *MAT_249 of LS-DYNA for a binder-stabilized UD-NCF by introducing three independent fiber families for the carbon rovings, glass fibers and stitching. They calibrated the model with tensile tests in the longitudinal and transverse direction as well as 45°-OATs. 30°- and 60°-OATs were used for the validation and a good agreement of the resulting forces and shear angles was achieved, but notable out-of-plane wrinkling was predicted. During application to hemispherical forming tests, the location and orientation of wrinkling were qualitatively predicted but resulting shear angles were overestimated and other strains were not further investigated in their study.

In the present work, a new hyperelastic macroscopic forming model for UD-NCF is proposed based on the preliminary work presented in [8], by introducing a strain energy density instead of non-linear stiffnesses and considering the previously omitted coupling of transverse tension and in-plane compression in the deformation behavior. The strain energy density is additively decomposed into individual components that are associated with different deformation modes, which is often assumed for engineering textiles [19,22–25,42]. For UD-NCF these are tension in fiber direction, transverse tension, shear deformation and in-plane compression perpendicular to the fiber rovings. A coupling between transverse tension and perpendicular compression is proposed to model large roving slippage if the rovings are under

<figure><img src="image_4.png" /><figcaption>Fig. 1. Unidirectional NCF UD300 from the (a) front and (b) back.</figcaption></figure>

compression. A new generalized approach with relatively few material parameters for the strain energy density and possible couplings between the deformation modes is used. The objective is to minimize the complexity of the model and thus facilitate its application to other materials. The model is parameterized using the forces and strains measured in different multiaxial deformation states obtained in OATs and quantitatively validated by forming tests on a hemisphere and tetrahedron punch shape.

## 2. Material

In this study, a unidirectional non-crimp fabric without binder is used, cf. Fig. 1. The fabric called UD300 is manufactured by Zoltek™ and produced from their PX35-50K continuous carbon fiber (CF) heavy tows, which are stitched together with a 76 dtex PES yarn in a Tricot pattern. The UD-NCF consists of a single layer of aligned CF tows with an areal weight of about 300 g m⁻² and thin glass fibers (GF) on the back as stitching basis and for improved handleability.

## 3. Hyperelastic membrane model for UD-NCF

### 3.1. Fundamentals of hyperelasticity adopted for UD-NCF

Hyperelastic approaches are defined by a total strain energy density function $W_{tot}$ as a potential for the resulting stress [19]

$$S = 2 \frac{\partial W_{tot}(C)}{\partial C}, \qquad (1)$$

where S is the second Piola–Kirchhoff stress tensor and C is the right Cauchy–Green deformation tensor. To account for anisotropic behavior, $W_{tot}$ can alternatively be expressed by a set of directional pseudo-invariants $I_i$ of the strain by

$$S = 2 \sum_{i}^{N} \frac{\partial W_{tot}}{\partial I_{i}} \frac{\partial I_{i}}{\partial C}. \qquad (2)$$

Thereby, physically interpretable pseudo-invariants are commonly chosen based on the principal orientations $a_0$ and $b_0$, cf. Fig. 2 a, to relate them to the observed deformation modes of the material. The directional pseudo-invariants [43]

$$I_4 = a_0 \cdot C \cdot a_0, \qquad (3)$$

$$I_6 = a_0 \cdot C \cdot b_0, \qquad (4)$$

$$I_8 = b_0 \cdot C \cdot b_0, \qquad (5)$$

$$I_{10} = \frac{\pi}{2} - \psi_{12} = \arccos(a_0 \cdot b_0) - \arccos\left(\frac{I_6}{\sqrt{I_4 I_8}}\right), \qquad (6)$$

$$I_{11} = \sqrt{I_8} \sin(\psi_{12}) = \sqrt{I_8} \sin\left(\arccos\left(\frac{I_6}{\sqrt{I_4 I_8}}\right)\right), \qquad (7)$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. (a) Schematic illustration of the orientation of pseudo-invariants in the initial and current configuration; (b) Generalized hyperelastic energy and its derivative; (c) Coupling functions with schematic illustrations of its parameters influence; (d) Schematic representation of different deformation modes (adapted from [41]) and the corresponding resulting pseudo-invariants.</figcaption></figure>

are introduced, where $I_4$ and $I_8$ are the quadratic stretches along $a_0$ and $b_0$ respectively, $I_6$ is the shear strain, $I_{10}$ is the shear angle $\gamma_{12}$, $\psi_{12}$ is the current angle between $a$ and $b$, and $I_{11}$ is the stretch in the derived direction $b^\perp$. The derived direction $b^\perp$ is the component of $b$ that remains perpendicular to the carbon fiber direction $a$ for arbitrary angles $\psi_{12}$, cf. Fig. 2 a. For UD-NCF the pseudo-invariants are related to tension in fiber direction ($I_4$), transverse tension ($I_8$), shear deformation ($I_6$) and in-plane compression perpendicular to the fiber rovings ($I_{11}$). The shear strain ($I_6$) is chosen instead of the shear angle ($I_{10}$), to ensure independence of the shear deformation from the stitching strains, cf. Eq. (6). Combining Eq. (2) with the invariants (Eq. (3)–(7)) and their partial derivatives A yields the expression

$$S = 2 \left( \frac{\partial W_{\text{tot}}}{\partial I_4} \frac{\partial I_4}{\partial C} + \frac{\partial W_{\text{tot}}}{\partial I_6} \frac{\partial I_6}{\partial C} + \frac{\partial W_{\text{tot}}}{\partial I_8} \frac{\partial I_8}{\partial C} + \frac{\partial W_{\text{tot}}}{\partial I_{11}} \frac{\partial I_{11}}{\partial C} \right) \quad (8)$$

$$= 2 \left( \frac{\partial W_{\text{tot}}}{\partial I_4} \frac{\partial I_4}{\partial C} + \frac{\partial W_{\text{tot}}}{\partial I_6} \frac{\partial I_6}{\partial C} + \frac{\partial W_{\text{tot}}}{\partial I_8} \frac{\partial I_8}{\partial C} \right) \quad (9)$$

$$+ 2 \left( \frac{\partial W_{\text{tot}}}{\partial I_{11}} \left( \frac{\partial I_{11}}{\partial I_4} \frac{\partial I_4}{\partial C} + \frac{\partial I_{11}}{\partial I_6} \frac{\partial I_6}{\partial C} + \frac{\partial I_{11}}{\partial I_8} \frac{\partial I_8}{\partial C} \right) \right)$$

$$= 2 \left( \frac{\partial W_{\text{tot}}}{\partial I_4} a_0 \otimes a_0 + \frac{\partial W_{\text{tot}}}{\partial I_8} b_0 \otimes b_0 \right) + \frac{\partial W_{\text{tot}}}{\partial I_6} (a_0 \otimes b_0 + b_0 \otimes a_0) \\ + 2 \frac{\partial W_{\text{tot}}}{\partial I_{11}} \frac{1}{\sqrt{I_4 I_8 - I_6^2}} \left( \frac{I_6^2}{2I_4^{3/2}} a_0 \otimes a_0 + \frac{-I_6}{\sqrt{I_4}} a_0 \otimes b_0 + \frac{\sqrt{I_4}}{2} b_0 \otimes b_0 \right). \quad (10)$$

The next step is the definition of the total strain energy density function $W_{\text{tot}}$ based on the observed deformation behavior of the UD-NCF.

## 3.2. General assumptions

UD-NCFs deform under large shear strains superimposed with transverse tensile deformation and in-plane compression perpendicular to the carbon fiber orientation. This behavior was experimentally observed during the characterization of the membrane behavior based on off-axis-tension tests (OATs) [6,40,41,44] as well as during forming tests [4,6,45,46]. Therefore, UD-NCF undergoes not only pure or simple shear, but also mixed deformation modes with different amplitudes and signs of multiaxial strains, cf. Fig. 2 d. The deformation of individual components is assumed to be directly related to the proposed invariants:

* Roving slippage parallel to the fiber direction results in shear $I_6$ as well as increased transverse tensile deformation $I_8$, as demonstrated during simple shear or the mixed-modes 2 & 3 in Fig. 2 d.

* Transverse tensile deformation $I_8 > I_8^0 = 1$ additionally indicates tension in the stitching pattern (initially aligned with $b_0$), assuming limited slippage between stitching and rovings.

* It is assumed that the stitching only bears tensile stresses $I_8 > I_8^0$, while in-plane compaction takes place via the rovings and is perpendicularly oriented to them $I_{11} < I_{11}^0 = 1$.

* The perpendicular invariant $I_{11}$ quantifies the distance between fiber rovings, indicating their compression $I_{11} < I_{11}^0$ as well as their widening or gap formation $I_{11} > I_{11}^0$.

In our preliminary study [8] we have demonstrated that a decoupled consideration of the deformation modes is not sufficient to model the
membrane behavior as experimentally demonstrated during OATs [41] (to be discussed further in Section 4). Significant roving slippage (high $I_6$&$I_8$) occurs in areas of high roving compression (mixed-mode 2, cf. Fig. 2 d), compared to cases where tensioning of the stitching limits the shear deformation (mixed-mode 3, cf. Fig. 2 d). A coupling between the transverse tensile ($I_8$) and perpendicular compressive behavior ($I_{11}$) is introduced to model this interaction.

Based on the above assumption, it is assumed that the total strain energy density $W_{tot}$ can be decomposed additively into tension of the rovings $W_4(I_4)$, shear deformation $W_6(I_6)$, transverse in-plane tension coupled to perpendicular compression $W_8(I_8, I_{11})$ and in-plane compaction perpendicular to the fiber rovings $W_{11}(I_{11})$ according to:

$$W_{tot}(C, a_0, b_0) = W_4(I_4) + W_6(I_6) + W_8(I_8, I_{11}) + W_{11}(I_{11}) \quad (11)$$

The constitutive equations for the individual energy components and coupling are formulated in a parameterized way in the following, to simplify the parameterization process.

### 3.3. Generalized strain energy density function

Piecewise functions [8,22] or polynomial functions [19,21,25] are usually chosen to define the strain energy density and subsequently stress–strain relation. Polynomial functions are highly adaptable and commonly fitted to analytical models that directly derive the stress–strain relation from experimental tests. Piecewise functions prescribe a characteristic curve, whereby the parameters have a localized impact on the trend and are thus often easier to interpret physically. Since the typical deformation modes of UD-NCF consist of superimposed strains, cf. Section 3.2, the respective contributions to the total energy cannot necessarily be identified individually during the characterization. Therefore, previous studies [6,8] have shown that an inverse parameterization process is required to simultaneously consider the force- and strain–displacement results of different characterization tests.

Based on the above assumptions, a generalized, piecewise function for the derivative of the strain energy density is chosen as basis for the individual contributions in Eq. (11)

$$\frac{\partial W_{i}^{G}(I_{i})}{\partial I_{i}} = \begin{cases} C_{i}^{\text{ini}} I_{i}^{a} & , [0 \le I_{i}^{a} < I_{i}^{t}] \\ S_{i}^{\text{ini}} + (C_{i}^{\text{mid}}(I_{i}^{a} - I_{i}^{t}) + \Delta S_{i}^{\text{mid}}) \times (1 - e^{-c_{i}^{t}(I_{i}^{a} - I_{i}^{t})}) & , [I_{i}^{t} \le I_{i}^{a} < I_{i}^{\text{lock}}] \\ S_{i}^{\text{ini}} + (C_{i}^{\text{mid}}(I_{i}^{a} - I_{i}^{t}) + \Delta S_{i}^{\text{mid}}) \times (1 - e^{-c_{i}^{t}(I_{i}^{a} - I_{i}^{t})}) + C_{i}^{\text{lock}}(I_{i}^{a} - I_{i}^{\text{lock}})^{2} & , [I_{i}^{\text{lock}} \le I_{i}^{a}] \end{cases} \quad (12)$$

with $I_i^a = |I_i - I_i^0|$, $I_i^t = \frac{S_{i}^{\text{ini}}}{c_{i}^{\text{ini}}}$ and $c_i^t = \frac{c_{i}^{\text{ini}}}{\Delta S_{i}^{\text{mid}}}$. The corresponding generalized strain energy density $W_i^G(I_i)$ is listed in Appendix B and the generalized curves are shown in Fig. 2 b. The formulation was chosen based on its high adaptability with six, easy to physically interpret parameters. It consists of an initial stress limit $S_{i}^{\text{ini}}$ and initial modulus $C_{i}^{\text{ini}}$. The modulus increases or decreases to a middle modulus $C_{i}^{\text{mid}}$ over a stress range $\Delta S_{i}^{\text{mid}}$ and increases again with $C_{i}^{\text{lock}}$ after $I_{i}^{\text{lock}}$. The transition invariant $I_i^t$ and exponent $c_i^t$ are determined to achieve an overall $C^1$-smooth function.

This generalized approach can represent the typical three-stage response ($C_{i}^{\text{mid}} < C_{i}^{\text{ini}}$) during shear or elongation [19,21,22] as well as the exponential character ($C_{i}^{\text{mid}} > C_{i}^{\text{ini}}$) typically obtained during compaction of fabrics [6,47]. The approach is based on the absolute invariant $I_i^a$ due to the exponential function, and the sign for tension and shear must be considered separately by

$$W_4(I_4) = \operatorname{sgn}(I_4 - I_4^0) \cdot W_4^G(I_4) \quad \text{and} \quad W_6(I_6) = \operatorname{sgn}(I_6 - I_6^0) \cdot W_6^G(I_6). \quad (13)$$

The strain energy contribution due to the perpendicular invariant $I_{11}$ only accounts for in-plane perpendicular compression

$$W_{11}(I_{11}) = \begin{cases} -W_{11}^{G}(I_{11}), & [I_{11} < I_{11}^{0}] \\ 0, & [I_{11}^{0} \le I_{11}] \end{cases} \quad (14)$$

The transverse strain energy density contribution $W_8(I_8, I_{11})$ is detailed separately in the following section due to the coupling to the perpendicular invariant $I_{11}$.

### 3.4. Transverse tension - in-plane compression coupling

In-plane compression in the perpendicular direction ($I_{11} < I_{11}^0$) facilitates large roving slippage (large $I_8$) [8,41]. A coupling function is introduced to model this behavior:

$$k_{11}(I_{11}) = \begin{cases} (1 - k_{11}^{\text{min}}) e^{-s_{11}(I_{11}^{a})^{c_{11}^{\text{mid}}}} + k_{11}^{\text{min}}, & [I_{11} \le I_{11}^{0}] \\ 1, & [I_{11}^{0} \le I_{11}] \end{cases} \quad (15)$$

and the impact of its parameters on the characteristic trend is shown in Fig. 2 c. The threshold slippage–compression-parameter $k_{11}^{\text{min}}$ is only allowed to attain values between 1 and 0, representing a maximum reduction of the strain energy density $W_8$ for large in-plane compression. The coefficient $s_{11}$ characterizes the intensity of the coupling and the exponent $e_{11}$ controls the onset of the coupling. The transverse strain energy density contribution is given by

$$W_8(I_8, I_{11}) = \begin{cases} 0, & [I_8 < I_8^0] \\ k_{11}(I_{11}) \cdot W_8^G(I_8), & [I_8^0 \le I_8] \end{cases} \quad (16)$$

A multiplicative coupling was chosen to take advantage of the generalized strain energy density function (Eq. (12)) and simplify the transferability to new materials or couplings in the future.

### 3.5. Implementation for forming simulations in ABAQUS/EXPLICIT

The following section details the implementation of the membrane behavior described above as well as the bending and interface behavior in the commercially available FE solver ABAQUS/EXPLICIT, which is applied for forming simulations in this work.

**Membrane behavior.** The proposed hyperelastic model is implemented by means of a user-defined material behavior (VUMAT) and applied to M3D3 membrane elements with their edges aligned in fiber direction to avoid numerical intra-ply locking [48,49]. Based on previous studies [6, 8], the material parameters are identified with an inverse approach by simultaneously considering the force- and strain–displacement results of three off-axis-tension tests (OATs). Details of the parameterization process for the membrane behavior are further described in Section 4.

**Bending behavior.** The bending behavior of UD-NCF needs to be considered in a decoupled fashion during forming due to the fabric's low transverse shear stiffness [6]. Since it is not the focus of this work, an existing approach for the bending behavior is used that was previously developed for thermoplastic tapes [50] and already successfully applied for UD-NCF [6] as well as woven fabrics [22]. Bending is modeled with a hypoelastic approach, that accounts for large shear deformation. It considers an orthotropic elastic bending stiffness in the resulting non-orthogonal material frame. The bending stiffness in roving direction $B_{11} = 4.06$ N mm and in transverse direction $B_{22} = 0.018$ N mm have been determined based on standard cantilever tests. The approach is implemented in a VUGENS subroutine for a user-defined integration scheme over the thickness of S3R shell elements, which are superposed to the membrane elements via shared nodes.

**Interface behavior.** Tool-ply and ply-ply contact is modeled via the built-in general contact algorithm in ABAQUS/EXPLICIT. The friction coefficients of $\mu^{\text{tool-ply}} = 0.147$ and $\mu^{\text{ply-ply}} = 0.336$ were determined with a sled-based pull-over setup [51] previously used for woven fabrics [52] according to ASTM standard D1894 [53].
<figure><img src="image_6.png" /><figcaption>Fig. 3. (a) Experimental setups and (b) numerical models for the parameterization with off-axis-tension tests (OATs).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison of the experimentally measured (left) and predicted (right) contour plots of the (a) 30°-, (b) 45°- and (c) 60°-OAT at a displacement of 80 mm.</figcaption></figure>

## 4. Parameterization of the membrane model with off-axis tension tests

Off-axis-tension tests (OATs) with three different bias angles (30°, 45° and 60°) are used to parameterize the membrane behavior [6,8].

**Experimental results.** The experimental results are detailed in Schäfer et al. [41] and the setups are shown in Fig. 3 a. A regular pattern of white dots with a distance of 10 mm was applied to the specimen for a DIC measurement of the macroscopic strains. Three distinct types of zones are identifiable independent of the bias angle: Rovings in zones “C” are clamped along one edge and barely deform, rovings in zones “A” are free at both ends resulting in the highest strains, and rovings in the secondary zones “B” are subjected to intermediate strains depending on their exact location. The OATs impose different ratios of large superimposed shear, transverse tension and perpendicular roving compression, making them suitable to determine the necessary material parameters for the hyperelastic model described in Section 3. In addition to the force-displacement curves, averaged strain-displacement curves over elements in the main deformation zone (MDZ) (highlighted in Fig. 3) are calculated.

**Observations during the parameterization process.** Parameter studies were performed, aiming at the best overall agreement of forces, local strains in the MDZ (cf. Fig. 6) as well as global strains (cf. Fig. 4) for the three OATs simultaneously. Therefore, equivalent strain measures to the pseudo-invariants were calculated by $E_{22} = I_8 - I_8^0$, $\gamma = I_{10}$ and $E_{\perp} = I_{11} - I_{11}^0$. The simulation models are shown in Fig. 3 b and each test is modeled by means of 8200 triangular elements with a

computation time of approx. 4 min on an Intel Core i5-8265U with 4 cores and 16 GB RAM. The tensile modulus $C_4^{\text{ini}}$ is assumed to be a sufficiently high constant value because no significant tensile strains were measured experimentally, as an alternative to a numerically complex inextensibility condition. The identification of the remaining material parameters is performed manually, by gradually introducing the parameters impacting the initial ($C_i^{\text{ini}}$, $S_i^{\text{ini}}$), middle ($\Delta S_i^{\text{mid}}$, $C_i^{\text{mid}}$) and locking sections ($I_i^{\text{lock}}$, $C_i^{\text{lock}}$) of the generalized approach as well as the coupling function ($k_{11}^{\text{min}}$, $s_{11}$, $e_{11}$), cf. Fig. 2 b and c. The stages of the parameterization process are not presented, but the most relevant findings are summarized in the following paragraph.

Each energy contribution influences the resulting forces and strains with varying intensity during different displacement ranges of the different OATs. In general, $W_{11}$ mainly influences the overall relative necking of the specimens and the compressive strain distribution, while only slightly impacting the forces. Transverse tension $W_8$ and shear $W_6$ on the other hand have a decisive influence on the force-displacement relation for all OATs, while establishing the order of magnitude for the strains depending on their relative ratio. The behavior during the 45°-OAT is mostly dominated by shear deformation, compared to the 30°- and 60°-OAT in which proportionally a higher amount of transverse tensile deformation occurs. Therefore, initial parameters for $W_6$ are determined based on the 45°-OAT and then $W_8$ is used to adapt the membrane model to the other orientations in the 30°- and 60°-OAT. The locking parts of $W_6$ and $W_8$ are not utilized in this work, since their introduction always led to a clear overestimation of the forces. However, the perpendicular compressive locking of $W_{11}$ is crucial to prevent a collapse of highly sheared elements and therefore
<figure><img src="image_8.png" /><figcaption>Fig. 5. Contributions to the strain energy density.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Overview of forces and strains in the MDZ of the off-axis-tension tests (OATs) with a bias direction of (a) 30°, (b) 45° and (c) 60°. Color bands indicate standard deviation.</figcaption></figure>

acts as an alternative to the shear locking mechanism for unidirectional fabrics. The identified material parameters are listed in Appendix D in Table D.1. The contributions to the total strain energy density due to shear $W_6$, perpendicular in-plane compression $W_{11}$ and transverse tension $W_8$ are shown in Fig. 5.

Final parameterized membrane behavior in the OATs. The general deformation behavior corresponds well between the simulations and experiments as evident from Fig. 4. The resulting forces and averaged shear angles over the MDZ are also in good agreement with the experimental results, cf. Fig. 6. During the 30°-OAT (cf. Fig. 6 a), the slope of the
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Schematics and dimensions of the utilized shapes: (a) hemisphere and (b) tetrahedron. (c) Position of local evaluation zones to calculate averaged strains over tool displacement.</figcaption></figure>

*Source:* adapted from [45].

forces in the simulation is too high at a displacement of 25 mm and flat-
tens out again from 30 mm because at this point the coupling strongly
reduces W<sub>8</sub>. The transverse tensile strains E<sub>22</sub> in the MDZ of the 30°-
OAT are overestimated for large displacements because the specimen's
main deformation concentrates in this zone and is homogeneously
distributed during simulation. In the experiments, the large transverse
tensile strains are more inhomogeneous and limited to slippage at the
transition between zones A and B2, cf. Fig. 4 a. Additionally, notable
shear and compressive strains also occur outside the MDZ due to the
lateral contraction of the stitching during the experiments. Both of
these effects cannot be modeled with a local macroscopic approach that
assumes homogeneous material properties over the whole specimen.

The perpendicular compressive strains are slightly underestimated
in the MDZ for all OATs. However, a further reduction of W<sub>11</sub> resulted
in an overestimated necking of the outer contours, so a compromise
of both objectives was chosen. Nevertheless, the simulation approach
was able to model the perpendicular compression in the MDZ of the
60°-OAT, which has not been achieved by previous modeling ap-
proaches [6,8]. In the experiments of the 45°- and 60°-OAT, positive
E<sub>22</sub>- and E<sub>⊥</sub>-strains occur along the outer edges of the specimens, cf.
Fig. 4. The stitching has slipped inward along the rovings in these
areas during the tests, resulting in a reduced transverse stiffness and
small gapping. This effect cannot be reproduced by the macroscopic
approach.

*Summary.* The proposed hyperelastic membrane approach achieves good agreement with the various ratios of superimposed shear, transverse tension and perpendicular compression strains imposed by the different OATs in the experiments. The deformation behavior is better predicted within the main deformation zone compared to the secondary deformation zones. This is mainly due to the complex behavior of the stitching during the experiments, as it can slip relative to the rovings and distribute loads over large areas due to lateral contraction. This behavior cannot be modeled with the proposed approach due to the limitations of a local macroscopic approach.

5. Application and validation

5.1. *Setup of forming tests*

The proposed hyperelastic membrane model is applied to forming
simulations of two double-curved geometries of a hemisphere and tetra-
hedron. Each layer is modeled by means of 14,400 triangular elements
with a computation time of approx. 20 min per layer on an Intel Core
i5-8265U with 4 cores and 16 GB RAM. The results are compared to
experimental forming tests performed at INSA Lyon, whose setup was
previously used for the validation of simulation approaches for woven
fabrics [18] or biaxial NCF [38]. The experimental results are detailed
in [45] and the schematic setups are shown in Fig. 7 a & b. Layers
of fabric are placed between an open blankholder and die made from
transparent PMMA material and deformed with a punch. Two single-
layer ([0°], [90°]) and a double-layer setup ([0°/90°]) were used for
the hemisphere test and three single-layer setups ([0°], [45°], [90°])
for the tetrahedron test. The forming depths are 75 mm and 95 mm
with blankholder weights of 1.7 kg and 6.0 kg for the hemisphere and
tetrahedron, respectively. A similar dot pattern to that used for the
OATs (cf. Fig. 3 a) and the same 2D DIC-based strain measurement
method (detailed in [41]) were applied to one half of the specimens.
The edges of the openings of the otherwise transparent lower dies
appear white on the images, with a similar contrast to the white dots
of the applied pattern. This results in a loss of correlation for points
passing over these edges, which becomes obvious in strongly deformed
elements, cf. Fig. 9. The 2D measurement has the disadvantage that
only the strains in the flat area can be measured accurately. In return,
the strains can be determined over the entire punch displacement.
Furthermore, there is no need for a specialized 3D measuring system,
which often requires a way to remove the textile from the mold and is
more suitable for binder-stabilized fabrics [7,17,31]. For a quantitative
evaluation of the strain development during the tests, the strains are
averaged over three local zones in areas of high deformation for each
shape, cf. Zones A, B and C in Fig. 7 c.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Outer contours (all in mm) of the (a) hemisphere and (b) tetrahedron forming tests.</figcaption></figure>

## 5.2. Hemisphere forming

The outer contours of the fabric layers for a displacement of $u = 75$ mm of the hemispherical punch are shown in Fig. 8 a. The contour is used to evaluate the overall deformation and often represents a first validation step for simulation approaches based on the material draw-in [17,37,38]. A strongly asymmetrical behavior unique to UD-NCF is observed, compared to the typically more symmetrical behavior previously reported for biaxial fabrics [18,38]. The results of the proposed approach are in good agreement with the experiments for the [0°]- and [45°]-test. The simulative contours of the double-layer tests nearly match the results of two superposed single-layer contours, because the individual layers deform mostly unconstrained due to the low blankholder weight. Similar behavior is observed during the experimental tests.

The experimental and simulative strains are shown in Fig. 9 for a more detailed comparison. The general orientation and location of areas with high shear deformation $\gamma_{12}$ are predicted well for all investigated setups. Considering the local shear angles in zones A, B and C, the magnitudes are also quantitatively in good agreement but with a tendency to be slightly underestimated, cf. Fig. 10 a. Comparing the [0°]-test and the 0°-layer in the [0°/90°]-test, the superposed 90°-layer reduces shear in zone A and increases it towards zones B and C both in simulation and experiment.

Local shear bands in the transverse direction (21-shear) are observed in Fig. 9 for all setups at the boundaries of the high shear area starting at the radius of the hemisphere, similar to those already reported for the model proposed by Schirmaier et al. [6]. In the experiments, shear deformation primarily occurs parallel to the carbon fiber rovings (12-shear), cf. Fig. 2 d, since the in-plane bending stiffness of the rovings prevents transverse in-plane shear. This distinction of the asymmetrical shear behavior cannot be achieved with the proposed macroscopic approach based on Cauchy mechanics. The introduction of an in-plane fiber bending stiffness could account for this phenomenon, but would require the calculation of higher strain gradients based on

generalized continua approaches [54,55] which are not yet realized for ABAQUS.

The negative transverse and perpendicular strains are well predicted within the areas of high shear deformation but overestimated in the outer areas. Especially in the [0°]-setup, positive strains in combination with small gaps between the rovings are observed along the lower edge during the experiments, which is similar to the observations during the 45°- and 60°-OATs and cannot be reproduced by the macroscopic simulation approach (cf. Section 4). The positive $E_{22}$- and $E_{\perp}$-strains in the center of the hemisphere during both single-layer tests are predicted well by the simulation. In the double-layer test, the additional superposed layer prevents positive strains in the center, which is captured by the model.

## 5.3. Tetrahedron forming

The outer contours of the fabric layers for a displacement of $u = 95$ mm of the tetrahedral punch are shown in Fig. 8 b. The general shapes of the simulative contours are in good agreement with the experiments. However, the material draw-in transverse to the initial fiber orientation is underestimated for all investigated setups. Additionally, large positive $E_{22}$- and $E_{\perp}$-strains are predicted in the tetrahedron's center during the simulation as opposed to the experimental measurement, cf. Fig. 11. This observation can be attributed to two reasons. First, the 2D-DIC method significantly underestimates the strains in the three-dimensional regions of the shape (demonstrated in Fig. C.12 in Appendix C). Second, the pressure due to the blankholder increases friction between the glass fibers, rovings and stitching, which results in a higher stiffness in the transverse direction [1]. This effect is not accounted for during the parameterization based on OATs and would also require considering the impact of the pressure on the transverse behavior along the glass-fiber direction, which is not possible with the current plane-stress approach. Both effects are more pronounced for the tetrahedron test compared to the hemisphere test, due to the stronger distortion in the DIC results and the higher blankholder weight.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison of experimental and simulative strains for [0°] and [45°] single-layer as well as [0°/90°] double-layer setups in the hemisphere forming tests for a punch displacement of *u* = 75 mm with a zoom on local shear bands due to transverse 21-shear (ts).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Shear angle in local zones of the (a) hemisphere and (b) tetrahedron forming test. Color bands indicate standard deviation.</figcaption></figure>

The general shear behavior and location of areas with high deformation are predicted well for the investigated setups, except for the area around zone B in the [0°]-setup, cf. Fig. 11. The high shear is

localized around the tetrahedron's corner and does not extend towards the edges due to the underestimated material draw-in in the transverse direction. This is also reflected in the quantitative comparison of the
<figure><img src="image_14.png" /><figcaption>Fig. 11. Comparison of experimental and simulative strains for [0°], [45°] and [90°] single-layer setups in the tetrahedron forming tests for a punch displacement of *u* = 95 mm.</figcaption></figure>

local strains, cf. Fig. 10 b. The local shear angles in all zones and configurations are in acceptable agreement between the simulation and experiment, except for the B in the [0°]-test. However, the onset of transverse shear bands is again noticeable in multiple areas, especially close to the transition from the flat parts to the tetrahedron.

While mixed-mode 1 ($I_8$, $I_{11} < 1$, cf. Fig. 2 d) was the main deformation mode in the hemisphere-test, mixed-mode 2 ($I_{11} < 1$, $I_8 > 1$) additionally occurs during tetrahedron forming near zone C in the [0°]-, zone B in the [90°]- and between zones B & C in the [45°]-setup. This deformation mode indicates significant roving slippage which is more pronounced in the simulation than in the experiments ($E_{22}^{sim} > E_{22}^{exp}$). However, it highlights the necessity to take into account the different deformation modes during characterization.

## 6. Conclusion

Unidirectional non-crimp fabrics deform not only under pure or simple shear, but under mixed deformation modes. These consist of different ratios of large shear superimposed with transverse tensile deformation due to roving slippage and stitching extension as well

as in-plane compression perpendicular to the carbon fibers. To describe this complex membrane behavior, a new hyperelastic approach is presented and implemented in the FE solverABAQUS/EXPLICIT. Suitable directional pseudo-invariants are introduced and related to the different deformation modes. Generalized formulations for the strain energy density contributions and coupling function are chosen to simplify the parameterization process and transferability to other materials. The transverse tensile energy is coupled to the perpendicular compression, to consider large roving slippage under stitching and roving compression.

The proposed hyperelastic approach is successfully parameterized based on the forces and different ratios of superimposed strains during off-axis tension tests with bias angles of 30°, 45° and 60°. A good overall agreement is achieved and the location as well as orientation of the different deformation zones observed in the experiments are well captured. It turned out, that a perpendicular compressive locking is essential to prevent a collapse of highly sheared elements, while a traditional shear locking and a transverse tensile locking are not necessary. Finally, the model is applied to forming of a hemisphere and tetrahedron geometry in different configurations and quantitatively compared to experimental results. The globally strongly asymmetrical
deformation behavior unique to UD-NCFs is well predicted during all tests. Additionally, the local shear angle in zones of high deformation was generally in good agreement between the simulation and experiments over the entire punch displacement with the tendency to be slightly underestimated in some places.

The results were especially good for the hemisphere compared to the tetrahedron with its geometrically challenging corners and higher blankholder weight. The modeling of the transverse behavior could be extended in the future to reflect the increasing adhesion between stitching, glass fibers and carbon rovings under normal pressure. Especially since this potentially impacts the behavior along the whole glass fiber direction, which was demonstrated experimentally by Schirmaier et al. [1]. This behavior requires the development of a non-local approach to increase the transverse tensile stiffness depending on the normal pressure of all elements along the glass fiber direction. Another limitation of the proposed macroscopic approach based on Cauchy mechanics is the symmetrical shear behavior, while for UD-NCF an asymmetrical shear behavior mainly parallel to the rovings is observed. This leads to transverse shear bands in the simulation which are prevented during experiments due to the in-plane bending stiffness of the rovings. The introduction of an in-plane fiber bending stiffness based on higher strain gradients could account for this phenomenon.

The approach and parameterization strategy were validated for a UD-NCF with a symmetrical tricot stitching pattern. A direct transferability of the method to UD-NCFs with comparable tricot [4,9,32,40] or other symmetrical stitching patterns [1,6] is expected due to similar observed membrane deformation modes. For quasi-unidirectional NCFs with a limited amount of 90°-fibers [35], it may be necessary to consider the implications of higher transverse tensile and perpendicular compressive stiffnesses, which will result in a reduced relevance of the coupling due to less observed roving slippage. For asymmetrical stitching patterns [33], it may be necessary to further distinguish between the transverse direction to account for roving slippage and the stitching direction. Therefore, a new contribution to the strain energy density could be introduced based on a directional pseudo-invariant in the asymmetrical stitching direction, in conjunction with the generalized stress function. Thus, the generalizability of the proposed deformation mechanisms for UD-NCFs with other architectures and the transferability of the method remains to be investigated.

### CRediT authorship contribution statement

**Bastian Schäfer:** Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. **Dominik Dörr:** Writing – review & editing, Software, Methodology. **Ruochen Zheng:** Writing – review & editing. **Naim Naouar:** Writing – review & editing, Supervision, Project administration. **Luise Kärger:** Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization.

### 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.

### Data availability

Data will be made available on request.

### Acknowledgments

The authors would like to thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) and the French National Research Agency (ANR) for funding the collaborative project "Composite forming simulation for non-crimp fabrics based on generalized

continuum approaches" (DFG: 431354059, ANR: ANR-19-CE06-0031), which the presented work is carried out for. This work is also part of the Heisenberg project "Digitalization of fiber-reinforced polymer processes for resource-efficient manufacturing of lightweight components", funded by the DFG (project no. 455807141)

## Appendix A. Derivatives of pseudo-invariants

In the following, the relevant derivatives of the utilized pseudo-invariants are listed:

$$ \frac{\partial I_4}{\partial C} = a_0 \otimes a_0 \qquad (A.1) $$

$$ \frac{\partial I_6}{\partial C} = \frac{1}{2} (a_0 \otimes b_0 + b_0 \otimes a_0) \qquad (A.2) $$

$$ \frac{\partial I_8}{\partial C} = b_0 \otimes b_0 \qquad (A.3) $$

$$ \frac{\partial I_{11}}{\partial C} = \frac{\partial I_{11}}{\partial I_4} \frac{\partial I_4}{\partial C} + \frac{\partial I_{11}}{\partial I_6} \frac{\partial I_6}{\partial C} + \frac{\partial I_{11}}{\partial I_8} \frac{\partial I_8}{\partial C} \qquad (A.4) $$

$$ \frac{\partial I_{11}}{\partial I_4} = \frac{I_6^2}{2I_4^{3/2} \sqrt{I_4 I_8 - I_6^2}} \qquad (A.5) $$

$$ \frac{\partial I_{11}}{\partial I_6} = \frac{-I_6}{\sqrt{I_4}} \frac{1}{\sqrt{I_4 I_8 - I_6^2}} \qquad (A.6) $$

$$ \frac{\partial I_{11}}{\partial I_8} = \frac{\sqrt{I_4}}{2} \frac{1}{\sqrt{I_4 I_8 - I_6^2}} \qquad (A.7) $$

## Appendix B. Generalized strain energy density function

Based on Eq. (12), the generalized strain energy density function is given by

$$ W_I^G(I_i) = \begin{cases} \frac{1}{2} C_i^{\text{ini}} I_i^a \qquad , [0 \le I_i^a < I_i^t] \\ S_i^{\text{ini}} I_i^a + \frac{1}{c_i^t C_i^{\text{ini}}} e^{-c_i^t (I_i^a - I_i^t)} \\ \times \left( C_i^{\text{ini}} (\Delta S_i^{\text{mid}} + C_i^{\text{mid}} I_i^a) + C_i^{\text{mid}} (\Delta S_i^{\text{mid}} - S_i^{\text{ini}}) \right) \\ + C_i^{\text{ini}} I_i^a \left( 2C_i^{\text{ini}} \Delta S_i^{\text{mid}} + C_i^{\text{ini}} C_i^{\text{mid}} I_i^a - 2C_i^{\text{mid}} S_i^{\text{ini}} \right) - W_1 \qquad , [I_i^t \le I_i^a < I_i^{\text{lock}}] \\ S_i^{\text{ini}} I_i^a + \frac{1}{c_i^t C_i^{\text{ini}}} e^{-c_i^t (I_i^a - I_i^t)} \\ \times \left( C_i^{\text{ini}} (\Delta S_i^{\text{mid}} + C_i^{\text{mid}} I_i^a) + C_i^{\text{mid}} (\Delta S_i^{\text{mid}} - S_i^{\text{ini}}) \right) \\ + C_i^{\text{ini}} I_i^a \left( 2C_i^{\text{ini}} \Delta S_i^{\text{mid}} + C_i^{\text{ini}} C_i^{\text{mid}} I_i^a - 2C_i^{\text{mid}} S_i^{\text{ini}} \right) - W_1 \qquad , [I_i^{\text{lock}} \le I_i^a] \\ + \frac{1}{3} C_i^{\text{lock}} (I_i^a - \Delta I_i^{\text{lock}})^3 \end{cases} \qquad (B.1) $$

$$ \text{with } W_1 = \frac{2(C_i^{\text{ini}} + C_i^{\text{mid}}) \Delta S_i^{\text{mid}}^2 + 2C_i^{\text{ini}} \Delta S_i^{\text{mid}} S_i^{\text{ini}} - C_i^{\text{mid}} S_i^{\text{ini}}^2}{2C_i^{\text{ini}}^2} + \frac{1}{2} \frac{S_i^{\text{ini}}^2}{C_i^{\text{ini}}} $$

## Appendix C. Distortion error of the 2D-DIC

To demonstrate the significant distortion in the tetrahedron test due to the 2D-DIC, the resulting node coordinates of the simulation were parallel projected onto the xy-plane. Equivalent 2D-strains were calculated with the DIC algorithm used during the experimental tests [41]. The results for the perpendicular strain $E_{\perp}^{2D}$ are shown in Fig. C.12.

## Appendix D. Material parameters

See Table D.1.
<figure><img src="image_15.png" /><figcaption>Fig. C.12. Exemplary conversion from 3D strains to 2D-equivalent strains by planar projection for the simulation of the tetrahedron forming tests.</figcaption></figure>

<table><caption>Table D.1
Material parameters of the hyperelastic membrane approach.</caption><thead><tr><th></th><th>C<sup>ini</sup><sub>i</sub><br>MPa</th><th>S<sup>ini</sup><sub>i</sub><br>MPa</th><th>C<sup>mid</sup><sub>i</sub><br>MPa</th><th>ΔS<sup>mid</sup><sub>i</sub><br>MPa</th><th>C<sup>lock</sup><sub>i</sub><br>MPa</th><th>I<sup>lock</sup><sub>i</sub><br>MPa</th></tr></thead><tbody><tr><td>W<sub>4</sub></td><td>1000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>W<sub>6</sub></td><td>0.14</td><td>0.018</td><td>0.01</td><td>0.005</td><td>-</td><td>-</td></tr><tr><td>W<sub>8</sub></td><td>1.1</td><td>0.025</td><td>0.03</td><td>0.02</td><td>-</td><td>-</td></tr><tr><td>W<sub>11</sub></td><td>0.01</td><td>0.0025</td><td>0.001</td><td>0.4</td><td>0.045</td><td>10.0</td></tr><tr><td></td><td>k<sup>min</sup><sub>11</sub></td><td>s<sub>11</sub></td><td>e<sub>11</sub></td><td></td><td></td><td></td></tr><tr><td>k<sub>11</sub></td><td>0.4</td><td>10.0</td><td>2.0</td><td></td><td></td><td></td></tr></tbody></table>

## References

[1] Schirmaier FJ, Weidenmann KA, Kärger L, Henning F. Characterisation of the draping behaviour of unidirectional non-crimp fabrics (UD-NCF). Composites A 2016;80:28–38. http://dx.doi.org/10.1016/j.compositesa.2015.10.004.

[2] Galkin S, Kunze E, Kärger L, Böhm R, Gude M. Experimental and numerical determination of the local fiber volume content of unidirectional non-crimp fabrics with forming effects. J Composites Sci 2019;3(1). http://dx.doi.org/10.3390/jcs3010019.

[3] Kunze E, Galkin S, Böhm R, Gude M, Kärger L. The impact of draping effects on the stiffness and failure behavior of unidirectional non-crimp fabric fiber reinforced composites. Materials (Basel, Switzerland) 2020;13(13). http://dx.doi.org/10.3390/ma13132959.

[4] Ghazimoradi M, Trejo EA, Carvelli V, Butcher C, Montesano J. Deformation characteristics and formability of a tricot-stitched carbon fiber unidirectional non-crimp fabric. Composites A 2021;82:106366. http://dx.doi.org/10.1016/j.compositesa.2021.106366.

[5] Azzouz R, Allaoui S, Moulart R. Composite preforming defects: a review and a classification. Int J Mater Form 2021;48(3):177. http://dx.doi.org/10.1007/s12289-021-01643-7.

[6] Schirmaier FJ, Dörr D, Henning F, Kärger L. A macroscopic approach to simulate the forming behaviour of stitched unidirectional non-crimp fabrics (UD-NCF). Composites A 2017;102:322–35. http://dx.doi.org/10.1016/j.compositesa.2017.08.009.

[7] Kärger L, Galkin S, Kunze E, Gude M, Schäfer B. Prediction of forming effects in UD-NCF by macroscopic forming simulation – capabilities and limitations. In: Proceedings of the 24th international conference on material forming. 2021. http://dx.doi.org/10.25518/esaform21.355.

[8] Schäfer B, Haas S, Boisse P, Kärger L. Investigation of the membrane behavior of UD-NCF in macroscopic forming simulations. Key Eng Mater 2022;926:1413–22. http://dx.doi.org/10.4028/p-2977b4.

[9] Ghazimoradi M, Montesano J. Macroscopic forming simulation for a unidirectional non-crimp fabric using an anisotropic hyperelastic material model. Appl Compos Mater 2023;55(3):171. http://dx.doi.org/10.1007/s10443-023-10158-0.

[10] Syerko E, Comas-Cardona S, Binetruy C. Models of mechanical properties/behavior of dry fibrous materials at various scales in bending and tension: A review. Composites A 2012;43(8):1365–88. http://dx.doi.org/10.1016/j.compositesa.2012.03.012.

[11] Syerko E, Comas-Cardona S, Binetruy C. Models for shear properties/behavior of dry fibrous materials at various scales: a review. Int J Mater Form 2015;8:1–23. http://dx.doi.org/10.1007/s12289-013-1138-7.

[12] Boisse P, Colmars J, Hamila N, Naouar N, Steer Q. Bending and wrinkling of composite fiber preforms and prepregs. A review and new developments in the draping simulations. Composites B 2018;141:234–49. http://dx.doi.org/10.1016/j.compositesb.2017.12.061.

[13] Bussetta P, Correia N. Numerical forming of continuous fibre reinforced composite material: A review. Composites A 2018;113:12–31. http://dx.doi.org/10.1016/j.compositesa.2018.07.010.

[14] Liang B, Boisse P. A review of numerical analyses and experimental characterization methods for forming of textile reinforcements. Chin J Aeronaut 2020. http://dx.doi.org/10.1016/j.cja.2020.09.027.

[15] Boisse P, Akkerman R, Carlone P, Kärger L, Lomov SV, Sherwood JA. Advances in composite forming through 25 years of ESAFORM. Int J Mater Form 2022;15(3):99. http://dx.doi.org/10.1007/s12289-022-01682-8.

[16] Peng X, Ding F. Validation of a non-orthogonal constitutive model for woven composite fabrics via hemispherical stamping simulation. Composites A 2011;42(4):400–7. http://dx.doi.org/10.1016/j.compositesa.2010.12.014.

[17] Khan MA, Mabrouki T, Vidal-Sallé E, Boisse P. Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark. J Mater Process Technol 2010;210(2):378–88. http://dx.doi.org/10.1016/j.jmatprotec.2009.09.027.

[18] Chen B, Colmars J, Naouar N, Boisse P. A hypoelastic stress resultant shell approach for simulations of textile composite reinforcement forming. Composites A 2021;149(4):106558. http://dx.doi.org/10.1016/j.compositesa.2021.106558.

[19] Aimène Y, Vidal-Sallé E, Hagège B, Sidoroff F, Boisse P. A hyperelastic approach for composite reinforcement large deformation analysis. J Compos Mater 2010;44(1):5–26. http://dx.doi.org/10.1177/0021998309345348.

[20] Charmetant A, Orliac JG, Vidal-Sallé E, Boisse P. Hyperelastic model for large deformation analyses of 3D interlock composite preforms. Compos Sci Technol 2012;72(12):1352–60. http://dx.doi.org/10.1016/j.compscitech.2012.05.006.
[21] Peng X, Guo Z, Du T, Yu W-R. A simple anisotropic hyperelastic constitutive model for textile fabrics with application to forming simulation. Composites B 2013;52:275-81. http://dx.doi.org/10.1016/j.compositesb.2013.04.014.

[22] Poppe C, Dörr D, Henning F, Kärger L. Experimental and numerical investigation of the shear behaviour of infiltrated woven fabrics. Composites A 2018;114:327-37. http://dx.doi.org/10.1016/j.compositesa.2018.08.018.

[23] Yao Y, Huang X, Peng X, Liu P, Youkun G. An anisotropic hyperelastic constitutive model for plain weave fabric considering biaxial tension coupling. Text Res J 2019;89(3):434-44. http://dx.doi.org/10.1177/0040517517748495.

[24] Yao Y, Peng X, Gong Y. Influence of tension-shear coupling on draping of plain weave fabrics. J Mater Sci 2019;54(8):6310-22. http://dx.doi.org/10.1007/s10853-019-03334-w.

[25] Schäfer F, Werner HO, Henning F, Kärger L. A hyperelastic material model considering biaxial coupling of tension-compression and shear for the forming simulation of woven fabrics. Composites A 2023;165(4):107323. http://dx.doi.org/10.1016/j.compositesa.2022.107323.

[26] Kärger L, Schäfer F, Werner HO. Modeling multiaxial stress states in forming simulation of woven fabrics. Mater Res Proc 2023;28:357-66. http://dx.doi.org/10.21741/9781644902479-39.

[27] Creech G, Pickett A. Meso-modelling of non-crimp fabric composites for coupled drape and failure analysis. J Mater Sci 2006;41(20):6725-36. http://dx.doi.org/10.1007/s10853-006-0213-6.

[28] Lomov S, Ivanov D, Verpoest I, Zako M, Kurashiki T, Nakai H, et al. Meso-FE modelling of textile composites: Road map, data flow and algorithms. Compos Sci Technol 2007;67(9):1870-91. http://dx.doi.org/10.1016/j.compscitech.2006.10.017.

[29] Sirtautas J, Pickett AK, Lépicier P. A mesoscopic model for coupled drape-infusion simulation of biaxial non-crimp fabric. Composites B 2013;47:48-57. http://dx.doi.org/10.1016/j.compositesb.2012.09.088.

[30] Khiem VN, Krieger H, Itskov M, Gries T, Stapleton SE. An averaging based hyperelastic modeling and experimental analysis of non-crimp fabrics. Int J Solids Struct 2016. http://dx.doi.org/10.1016/j.ijsolstr.2016.12.018.

[31] Mallach A, Härtel F, Heieck F, Fuhr J-P, Middendorf P, Gude M. Experimental comparison of a macroscopic draping simulation for dry non-crimp fabric preforming on a complex geometry by means of optical measurement. J Compos Mater 2017;51(16):2363-75. http://dx.doi.org/10.1177/0021998316670477.

[32] Böhler P, Härtel F, Middendorf P. Identification of forming limits for unidirectional carbon textiles in reality and mesoscopic simulation. Key Eng Mater 2013;554-557:423-32. http://dx.doi.org/10.4028/www.scientific.net/KEM.554-557.423.

[33] Senner T, Kreissl S, Merklein M, Meinhardt J, Lipp A. A modular modeling approach for describing the in-plane forming behavior of unidirectional non-crimp-fabrics. Prod Eng 2014;8(5):635-43. http://dx.doi.org/10.1007/s11740-014-0561-z.

[34] Kärger L, Bernath A, Fritz F, Galkin S, Magagnato D, Oeckerath A, et al. Development and validation of a CAE chain for unidirectional fibre reinforced composite components. Compos Struct 2015;132:350-8. http://dx.doi.org/10.1016/j.compstruct.2015.05.047.

[35] Krogh C, Kepler JA, Jakobsen J. Pure and simple: investigating the in-plane shear kinematics of a quasi-unidirectional glass fiber non-crimp fabric using the bias-extension test. Int J Mater Form 2021;19(3-4):1483-95. http://dx.doi.org/10.1007/s12289-021-01642-8.

[36] Yu W-R, Harrison P, Long A. Finite element forming simulation for non-crimp fabrics using a non-orthogonal constitutive equation. Composites A 2005;36(8):1079-93. http://dx.doi.org/10.1016/j.compositesa.2005.01.007.

[37] Chen S, McGregor O, Harper LT, Endruweit A, Warrior NA. Defect formation during preforming of a bi-axial non-crimp fabric with a pillar stitch pattern. Composites A 2016;91:156-67. http://dx.doi.org/10.1016/j.compositesa.2016.09.016.

[38] Guzman-Maldonado E, Bel S, Bloom D, Fideu P, Boisse P. Experimental and numerical analyses of the mechanical behavior during draping of non-orthogonal bi-axial non-crimp fabric composite reinforcements. Mater Des 2022;218(3):110682. http://dx.doi.org/10.1016/j.matdes.2022.110682.

[39] Bel S, Hamila N, Boisse P, Dumont F. Finite element model for NCF composite reinforcement preforming: Importance of inter-ply sliding. Composites A 2012;43(12):2269-77. http://dx.doi.org/10.1016/j.compositesa.2012.08.005.

[40] Ghazimoradi M, Trejo EA, Butcher C, Montesano J. Characterizing the macroscopic response and local deformation mechanisms of a unidirectional non-crimp fabric. Composites A 2022;156(9):106857. http://dx.doi.org/10.1016/j.compositesa.2022.106857.

[41] Schäfer B, Zheng R, Naouar N, Kärger L. Membrane behavior of uni- and bidirectional non-crimp fabrics in off-axis-tension tests. Int J Mater Form 2023;16(6):68. http://dx.doi.org/10.1007/s12289-023-01792-x.

[42] Boisse P, Bai R, Colmars J, Hamila N, Liang B, Madeo A. The need to use generalized continuum mechanics to model 3D textile composite forming. Appl Compos Mater 2018;47(11). http://dx.doi.org/10.1007/s10443-018-9719-8.

[43] Holzapfel G. Nonlinear solid mechanics: A continuum approach for engineering. West Sussex, England: John Wiley & Sons Lt.; 2000.

[44] Trejo EA, Ghazimoradi M, Butcher C, Montesano J. Assessing strain fields in unbalanced unidirectional non-crimp fabrics. Composites A 2020;130:105758. http://dx.doi.org/10.1016/j.compositesa.2019.105758.

[45] Schäfer B, Zheng R, Colmars J, Platzer A, Naouar N, Boisse P, et al. Experimental analysis of the forming behavior of uni- and bidirectional non-crimp fabrics for different geometries. Composites B 2024. (under revision).

[46] Kunze E, Schwarz B, Weber T, Müller M, Böhm R, Gude M. Forming analysis of internal plies of multi-layer unidirectional textile preforms using projectional radiography. Procedia Manuf 2020;47:17-23. http://dx.doi.org/10.1016/j.promfg.2020.04.110.

[47] Schäfer B, Zheng R, Boisse P, Kärger L. Investigation of the compaction behavior of uni- and bidirectional non-crimp fabrics. Mater Res Proc 2023;28:331-8. http://dx.doi.org/10.21741/9781644902479-36.

[48] ten Thije R, Akkerman R. Solutions to intra-ply shear locking in finite element analyses of fibre reinforced materials. Composites A 2008;39(7):1167-76. http://dx.doi.org/10.1016/j.compositesa.2008.03.014.

[49] Hamila N, Boisse P. Locking in simulation of composite reinforcement deformations. analysis and treatment. Composites A 2013;53:109-17. http://dx.doi.org/10.1016/j.compositesa.2013.06.001.

[50] Dörr D, Schirmaier FJ, Henning F, Kärger L. A viscoelastic approach for modeling bending behavior in finite element forming simulation of continuously fiber reinforced composites. Composites A 2017;94:113-23. http://dx.doi.org/10.1016/j.compositesa.2016.11.027.

[51] Schäfer B, Naouar N, Kärger L. Investigation of the friction behavior of uni- and bidirectional non-crimp fabrics. Mater Res Proc 2024;41:540-8. http://dx.doi.org/10.21741/9781644903131-60.

[52] Poppe C, Dörr D, Kraus F, Kärger L. Experimental and numerical investigation of the contact behavior during FE forming simulation of continuously reinforced composites in wet compression molding. In: AIP conference proceedings, AIP Publishing; 2019, 020002. http://dx.doi.org/10.1063/1.5112507.

[53] ASTM Standard. Test method for static and kinetic coefficients of friction of plastic film and sheeting: D1894-14. 2014, http://dx.doi.org/10.1520/D1894-14.

[54] Boisse P, Hamila N, Madeo A. The difficulties in modeling the mechanical behavior of textile composite reinforcements with standard continuum mechanics of Cauchy. Some possible remedies. Int J Solids Struct 2018;154:55-65. http://dx.doi.org/10.1016/j.ijsolstr.2016.12.019.

[55] Steer Q, Colmars J, Naouar N, Boisse P. Modeling and analysis of in-plane bending in fibrous reinforcements with rotation-free shell finite elements. Int J Solids Struct 2021;222-223. http://dx.doi.org/10.1016/j.ijsolstr.2021.03.001.