# A solid-beam approach for mesoscopic analysis of textile reinforcements

Baptiste Lacroix, Julien Colmars, Auriane Platzer, Emmanuelle Vidal-Sallé

► To cite this version:

Baptiste Lacroix, Julien Colmars, Auriane Platzer, Emmanuelle Vidal-Sallé. A solid-beam approach for mesoscopic analysis of textile reinforcements. Composite Structures, 2025, 379, pp.119943. 10.1016/j.compstruct.2025.119943. hal-05424052

HAL Id: hal-05424052

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

Submitted on 18 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# A solid-beam approach for mesoscopic analysis of textile reinforcements

Baptiste Lacroix, Julien Colmars<sup>ID*</sup>, Auriane Platzer<sup>ID</sup>, Emmanuelle Vidal-Sallé<sup>ID</sup>

INSA Lyon, CNRS, LaMCoS, UMR5259, 69621 Villeurbanne, France

## ARTICLE INFO

**Keywords:**

Fibrous solid-beam
Mesoscopic analysis
Transverse isotropic hyperelasticity
Composite reinforcements
Textile modeling

## ABSTRACT

In composite design and manufacturing process, numerical mesoscopic analysis aims to offer accurate tools to predict fiber volume fraction, gaps, fiber disorientation and other defect or geometrical characteristic that could influence mechanical properties or permeability of the fibrous media.

This work proposes a new approach for the mesoscopic analysis of textiles, based on a fibrous solid-beam finite element and a transversely isotropic hyperelastic constitutive model adapted from previous work. The element is designed to simulate fiber yarns in textiles, and is based on the superposition of a solid element and virtual beams that take into account the bending behavior of fibrous materials. The element's potential is illustrated through numerical simulations on a single yarn, and the role of the material parameters is discussed in detail.

## 1. Introduction

### 1.1. Numerical simulation of fibrous reinforcements

The manufacturing process of composite parts with complex geometry (e.g. RTM [1]) require a good understanding of fibrous reinforcements behavior, in particular during forming step. The process raises questions such as the final fiber orientation after forming, the feasibility of a given geometry, or the occurrence of local defects. The issue at stake here is often the durability of the composite part in service. The quality of the composite part depends in part on our ability to control the manufacturing stages, particularly forming steps, and more specifically fiber density, which controls permeability during the injection stage [2–4].

Numerical simulation is one way of improving our understanding of the mechanisms at work in the forming process. Significant progress has been made in this field in recent years, in particular at mesoscopic and microscopic levels [5]. At mesoscopic level, yarn (fiber bundle) is supposed to be a continuous medium, whereas microscopic approaches tend to consider all the fibers present in textiles (i.e. around ten thousand per yarn), or virtual fibers (bundle of hundreds fibers) [6,7]. The mesoscopic scale now seems to be a good numerical compromise for modeling textile behavior: it allows modeling of yarn slippage and misalignment [8], yarn buckling and distortion, gaps between yarns [9], and allows multi-scale approaches, for instance in the objective of determining macroscopic behavior with meso-models [10,11].

### 1.2. A variety of mesoscopic modeling approaches

The recent review of [5] demonstrated that a large variety of models has been dedicated to fibrous reinforcements modeling at meso-scale. This subject raises issues relating to: the geometric description of the medium, whether or not based on experimental observations (often using micro-tomography [12–15]); the mechanical characterization of the continuous medium under consideration; and the numerical resolution methods used to simulate the forming process. This paper focus on the two last aspects: the mechanical behavior and the simulation methods at mesoscopic scale.

On the aspect of mesoscopic-scale behavior, the work presented here is a continuation of [16], itself inspired by the work of [17]. Although Charmetant's hyperelastic law has been widely used in recent work, its formulation and the determination of the associated material parameters still raise many questions, which we will explore in this article.

As far as numerical simulation tools are concerned, numerous finite element approaches have been proposed in recent years. They cover a wide range of methods, including the use of solid elements [18], shell elements [19] and beam elements [8,20]. Most of these models use first gradient mechanics, but some work have explored the interest of second gradient mechanics [21] which allows the material's bending stiffness to be taken into account, while decoupling it from other properties, which is typical of fibrous media [22].

To combine the advantages of the solid and beam approaches, some approaches explore beam-solid models (beam with cross-sectional

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: julien.colmars@insa-lyon.fr (J. Colmars).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Representation of the fibrous solid-beam finite element in its reference configuration.</figcaption></figure>

deformation) using enriched kinematics [23,24]. In this paper, we propose an approach based on solid-beam elements. This type of finite element has undergone numerous developments in recent decades; some of these strategies are detailed in the next section.

## 1.3. Solid-beam elements in literature

The scientific literature on solid-beam finite elements is fairly dense. Here, we focus on their applicability to fibrous media in the broadest sense, i.e. the modeling of slender fibrous structures (cables, wires, helical strands, etc.).

In this work, the term solid-beam refers to: elements based on a solid element geometric support, which must allow deformation of the cross-section transverse to the fiber direction, while using beam kinematics allowing bending stiffness to be considered. Solid-beam elements generally use refined kinematics [25–27], and sometimes iso-geometric analysis [21].

The idea of using a solid-beam approach is consistent with recent works using solid-shell for textile forming [28]. These type of elements are subjected to locking problems that are often solved using reduced integration of assumed-strain methods [29]. It was applied to both solid-shell elements [30] or solid-beam elements [26].

Very few solid-beam elements have been adapted to the context of fibrous media: some were developed in the assumption of small strain [31,32], and to our knowledge, the only approach including large strain is the one from [21]. Their approach can be considered to be specific to fibrous materials, since it uses an anisotropic hyperelastic behavior with a unique fiber direction, and because it also includes bending stiffness using second gradient approach. Although developed in the context of shells, the approach of [33] using Cosserat mechanics should be mentioned as an alternative strategy to take into account bending in fibrous media.

In order to overcome the difficulties inherent in second gradient approaches, and to facilitate contact management in existing codes, we propose to explore the possibility of a fibrous solid-beam element based on a rotation-free approach: the neighbor element method, which has been successfully used to reproduce in-plane bending in shells [34] will then be used as a substitute to a second-gradient formulation, in order to build a so-called fibrous solid-beam element.

## 1.4. Scope of the paper

In this paper, we propose a new fibrous solid-beam element: it is based on a solid element with a hyperelastic behavior, enriched with beams acting as fictive fibers. This element is designed to represent fiber yarns, defined as bundles of fibers aligned along a preferred direction. The development of the element and the constitutive model are presented in Section 2. Section 3 is devoted to the study of the influence of material parameters on roving behavior at the mesoscopic scale. Property ratios will be introduced to define different classes of

fiber bundle kinematic response. Finally, in Section 4 we propose a synthesis designed to help simulators and experimentalists in the choice of material parameters suitable for this type of medium.

## 2. Fibrous solid-beam finite element

In this study, yarns are meshed with only one element in their thickness (i.e. the smallest dimension). In textile models, yarns cross-section can have complex geometry (circular, ellipsoidal, lenticular, etc.). This paper focuses on a hexahedral solid-beam element; this implies that the cross-section of the yarn is represented as a quadrilateral when using a single element, but the geometry could be more precisely described by re-meshing in transverse direction. However, to concentrate on the material behavior and to demonstrate the capability of the fibrous solid-beam element itself, the following numerical examples will only involve one element in the cross-section.

### 2.1. Element geometry

The solid-beam element proposed here is composed of an 8-nodes isoparametric linear hexahedral element which supports four fictive 2-nodes fibers aligned with the preferential direction of the yarn. The fibrous element is represented in its referential state in Fig. 1. The four fictive fibers at the edges of the element represent a continuous distribution of fictive fibers across the thickness of the element. The assumption of the existence of fictive fibers is directly inherited from the specific fibrous element formulation developed in [35,36]. It was originally motivated by the inability of classical continuum mechanics to properly model fibrous media behavior, in particular the observed decoupling of tension and bending [22].

We can therefore consider the solid-beam fibrous element as the superposition of one hexahedral element and four fibers, which is consistent with the “modal” energy decomposition detailed hereafter. In addition, this approach allows for the decoupling of first and second gradient physics (bending), without using higher order elements.

### 2.2. Yarn mechanical behavior at mesoscopic scale

Previous research has shown, by observing the fiber volume fraction and the yarn’s cross-section shape on micro-tomography images [37], that the yarn can reasonably be considered as a transversely isotropic medium. Following these observations, a transversely isotropic hyper-elastic model was proposed by [16] to model yarn behavior in solid elements with a set of physical invariants. The same invariants will be used herein for the solid-beam element for compaction, distortion and longitudinal shear of the yarn, adapting the hyperelastic potentials to our needs.

Charmetant’s model follows Criscione’s assumptions and framework [17], in particular to calculate the constitutive basis ($M, N_1, N_2$) of any transformation, as defined by the deformation gradient tensor $F$. The constitutive basis is represented in Fig. 2, where $M$ is aligned with fiber direction and ($N_1, N_2$) are two vectors of the orthogonal plane. The latter are defined as the two orthogonal vectors which remain orthogonal under the deformation $F$. In practice, $N_1$ and $N_2$ are the eigenvectors of the inverse of the right Cauchy–Green tensor $C^{-1}$, reduced to its in-plane components, orthogonal to $M$. For the rest of this document, ($N_1, N_2$) will be referred to as the transverse plane. In the solid-beam element, the fiber direction $M$ is defined as the mean direction given by the four fictive fibers:

$$ \underline{M} = \frac{1}{4} \sum_{b=1}^{4} e_{f_b} \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Decomposition of transformation for a fibrous yarn.</figcaption></figure>

where $e_{f,b}$ is the direction of fiber b in the initial state, calculated from the initial position $\underline{X}_i$ of each node as (see Fig. 1):

$$ 
\begin{aligned}
e_{f1} &= \frac{\underline{X}_2 - \underline{X}_1}{\|\underline{X}_2 - \underline{X}_1\|} \\
e_{f2} &= \frac{\underline{X}_3 - \underline{X}_4}{\|\underline{X}_3 - \underline{X}_4\|} \\
e_{f3} &= \frac{\underline{X}_6 - \underline{X}_5}{\|\underline{X}_6 - \underline{X}_5\|} \\
e_{f4} &= \frac{\underline{X}_7 - \underline{X}_8}{\|\underline{X}_7 - \underline{X}_8\|}
\end{aligned} 
\quad (2) $$

The orthonormal direct triad $(\underline{M}, \underline{N}_1, \underline{N}_2)$ in the initial state is then transformed into the triad $(\underline{m}, \underline{n}_1, \underline{n}_2)$ in the current state by a rigid body rotation $\underline{Q}$:

$$ \underline{m} = \underline{Q} \cdot \underline{M} \quad (3a) $$

$$ \underline{n}_1 = \underline{Q} \cdot \underline{N}_1 \quad (3b) $$

$$ \underline{n}_2 = \underline{Q} \cdot \underline{N}_2 \quad (3c) $$

In addition, the constitutive basis is oriented thanks to the following criterion such that $\underline{N}_1$ (respectively $\underline{N}_2$) can be interpreted as the direction in which fiber density becomes the lowest (respectively the highest) in the current state:

$$ \underline{N}_1 \cdot \underline{C}^{-1} \cdot \underline{N}_1 \le \underline{N}_2 \cdot \underline{C}^{-1} \cdot \underline{N}_2 \quad (4) $$

To illustrate this, we represent in Fig. 3 an arbitrary deformation of a yarn composed of one compression and one rigid body rotation. The orientation of the constitutive basis, in the initial state, is aligned with the compression direction $\underline{N}_2$. The constitutive basis in the current state is deduced from Eq. (3); $\underline{n}_2$ (respectively $\underline{n}_1$) then gives the direction in which fiber volume fraction is indeed the highest (respectively the lowest).

Consequently, [17] introduced a specific decomposition of the deformation gradient $\underline{F} = \underline{Q} \cdot \underline{f}$ which reveals a so-called intrinsic deformation gradient tensor $\underline{f}$, removing the rigid body rotation $\underline{Q}$. This decomposition is also illustrated in Fig. 2.

Under these assumptions, and further developments listed in Appendix A, the intrinsic deformation gradient tensor $\underline{f}$ is expressed in the constitutive basis as follows:

$$ \underline{f} = \begin{bmatrix} \underline{f}_M & \underline{f}_{M1} & \underline{f}_{M2} \\ 0 & \underline{f}_{11} & 0 \\ 0 & 0 & \underline{f}_{22} \end{bmatrix}_{(\underline{M}, \underline{N}_1, \underline{N}_2)} \quad (5) $$

where the $\underline{f}_{ij}$ components will be of further use in the definition of modal invariants.

Following Charmetant, the strain energy density is decomposed into different modes as in Eq. (6). The four modes proposed by Charmetant

are conserved (tension $W_{ten}$, compaction $W_{comp}$, distortion $W_{dist}$ and shearing $W_{shear}$) and a specific bending energy $W_{bend}$ is added to build the solid-beam element:

$$ W = W_{bend} + W_{ten} + W_{comp} + W_{dist} + W_{shear} \quad (6) $$

As elaborated in the subsequent sections, the tensile and bending behaviors are encoded in the fictive fibers while the other deformation modes that affect the shape and volume of the yarn are encoded in the hexahedral element.

## 2.3. Modes in fictive fibers

### 2.3.1. Bending

Fictive fibers are modeled as Bernoulli beams and their curvature calculated thanks to a non-local method called Neighbor Element Method (NEM). NEM is based on relative rotations between the element of interest and its two neighbors. It only uses displacement degrees of freedom, which readily enforces the interaction between the fictive fibers and the hexahedral element. Historically, NEM was developed for stamping simulation, and more precisely for bending calculation in shells [38–40]. Since then, the method has been adapted to reproduce in-plane bending effects in shells [34,41], imitating the effect of second gradient.

The NEM method has been adapted herein to a 3D beam formulation. In three dimensions, the curvature is now treated as a two-components vector, whose components are expressed in an orthonormal basis associated to the local fictive fiber orientation. This basis is arbitrarily defined in the initial state as an orthonormal triad oriented by the local fictive fiber direction $e_{f,b}$ (see Eq. (2)). The two planes containing the local fictive fiber direction are the one in which NEM will be applied. To perform the projection in these planes, the local basis is tracked at each calculation step in order to determine the rotation between the current and initial state. In Fig. 4 an isolated fiber following an helical trajectory is represented. It is meshed with three elements (one element of interest and its two neighbors). An arbitrary orthonormal basis $(\underline{f}_{-10}, \underline{f}_{-20}, \underline{f}_{-30})$ is attached to the element of interest in its initial state, drawn in black dashed line, while another basis $(\underline{f}_1, \underline{f}_2, \underline{f}_3)$ is attached to its current state, drawn in solid black line. The two dashed blue and green lines are respectively the projection of the fiber (each element) in $(\underline{f}_{-10}, \underline{f}_{-20})$ and $(\underline{f}_{-30}, \underline{f}_{-10})$ planes. The standard NEM is then carried out in each of these planes, thus yielding the two components of the curvature vector $\underline{\chi}_{f,b}$.

This 3D NEM version was validated through comparative tests with Euler–Bernoulli beam element (B33 in Abaqus), using general beam section option to allow decoupling between tension and bending stiffness. A cantilever test was calculated on a single fiber initially aligned with x global referential direction. Displacement of the free tip is controlled with respectively 30 mm and 15 mm deflection in z and y global referential directions. Comparing the curvatures calculated by NEM and Abaqus, the discrepancy between the two simulations never exceeded 2% for large deformation.
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Initial state of a fiber bundle with its constitutive basis and completely random fiber distribution; (b) Current state of compressed yarn’s cross-section with its constitutive basis and particular fiber distribution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Bending calculation in 3D on a single fiber using Neighboring Element Method (NEM). Representation of initial (black dashed line) and current state (black solid line) basis used for NEM, with the two projections used for curvature components calculation.</figcaption></figure>

In the following, bending behavior is considered as elastic, and isotropic. The bending stiffness of a yarn is evenly distributed on each fictive fibers; in the case of one element in the cross section, it is then divided by four. By denoting $\chi_{fb}$ the curvature vector of a fictive fiber b in its current configuration and $K_{bend}$ the bending stiffness of the fictive fibers, we can express bending energy density as:

$$W_{bend} = \sum_{f=1}^{4} \frac{K_{bend}}{2} \|\chi_{fb}\|^2 \qquad (7)$$

### 2.3.2. Tension

As illustrated in Fig. 7, the tension of the yarn is supported by truss elements attached to the edges of the hexahedral element, which are aligned with the mean fiber direction. Since the fibers are quasi-inextensible, this mode is usually much stiffer than the others, in order to ensure an accurate kinematic response, representative of the experimental results. A homogeneous fiber extension is presented on Fig. 5, where $\lambda_{f1}^0 = 1$ (respectively $\lambda_{f1}$) is the elongation of the bottom right fiber in its initial (respectively current) state.

The bottom right fictive fiber connecting nodes 1 and 2 is aligned with $\underline{e}_{f1}$, a unitary vector calculated in Eq. (2) from $\underline{X}_1$ and $\underline{X}_2$, the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Elongation mode.</figcaption></figure>

position of the nodes in the initial state. Elongation $\lambda_{f1}$ is calculated thanks to the projection of the right Cauchy-Green tensor $\underline{C} = \underline{F}^T \underline{F}$ on the initial state fictive fiber direction $\underline{e}_{f1}$ as in Eq. (8). The same procedure is applied to every fibers, the only difference being the coordinate of fictive fiber tip nodes and the point where tensor $\underline{C}$ is calculated.

$$\lambda_{f1} = \sqrt{\underline{e}_{f1} \cdot \underline{C} \cdot \underline{e}_{f1}} \qquad (8)$$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Deformations modes encoded in the hexahedral element.</figcaption></figure>

The motivation for encoding tension in the truss elements instead of the solid element as originally proposed by Charmetant is to conform with the fibrous element formulation: the behavior is written at the fiber level, instead of at the yarn level. Therefore, at this scale, we can opt for a classical continuum elastic media modeling. Thus, the tension of the fictive fiber simply follows a neo-Hookean hyperelastic model [42], whose energy potential depends on the tensile stiffness of the fiber $C_{10}^f$ and the first invariant of the right Cauchy-Green tensor $I_1 = \text{tr}(\underline{C})$:

$$W_{\text{NH}} = C_{10}^f (I_1 - 3) \qquad (9)$$

Considering each fictive fiber extension as an incompressible uniaxial tensile problem of a rod, for which $I_1$ expression is easily deduced, we sum each contribution to write the total strain energy density at the yarn level as:

$$W_{\text{ten}} = \sum_{b=1}^{4} C_{10}^{f} (2\lambda f_{b}^{2} + \lambda f_{b}^{-1} - 3) \qquad (10)$$

The stiffness $C_{10}^f$ in Eq. (10) can take any value but in fibrous media it is often much larger than the other stiffness components. In particular for carbon or glass fibers, fibers are considered inextensible in many simulations. Its value is then of little experimental significance; rather, it is a numerical value whose order of magnitude can be adjusted. With all these assumptions, the energy density associated to the tension mode is calculated by

$$W_{\text{ten}} = C_{10} \sum_{b=1}^{4} 2\lambda f_{b}^{2} + \lambda f_{b}^{-1} - 3. \qquad (11)$$

## 2.4. Modes associated to the hexahedral element

This section develops the last three modes considered for the yarn's behavior, namely compaction, distortion and longitudinal shearing, which are encoded in the hexahedral element, as represented in Fig. 6. The definition of the modes is built upon Charmetant's decomposition of the intrinsic deformation gradient tensor $\underline{f}$ [16]:

$$\underline{f} = \underline{f}_{\text{elong}} \cdot \underline{f}_{\text{comp}} \cdot \underline{f}_{\text{dist}} \cdot \underline{f}_{\text{shear}} \qquad (12)$$

In our case, the elongation mode defined by $\underline{f}_{\text{elong}}$ has been substituted by the fictive fiber's elongation presented in Section 2.3.2. Such consideration is made possible thanks to modal independence in Charmetant's law formulation [16]. In the following, each remaining modal deformation gradient tensor $\underline{f}$ will be expressed with the components in the constitutive basis introduced in Eq. (5).

### 2.4.1. Compaction

Compaction mode, represented in Fig. 6(a), is associated to a change in the cross-sectional area and can be interpreted as the variation of the fiber volume fraction of the yarn.

In this regard, it can be expressed as the spherical part of the strain tensor reduced to its two-dimensional form in the transverse plane

($\underline{N}_1$, $\underline{N}_2$). With the intrinsic deformation gradient tensor expression in Eq. (5), the compaction part can be written as

$$\underline{f}_{\text{comp}} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & f_{11}f_{22} & 0 \\ 0 & 0 & f_{11}f_{22} \end{bmatrix}_{(\underline{M}, \underline{N}_1, \underline{N}_2)} \qquad (13)$$

Its invariant

$$I_{\text{comp}} = \frac{1}{2} \ln (f_{11}f_{22}) \qquad (14)$$

quantifies the cross-sectional area ratio between initial and current state, components $f_{11}$ and $f_{22}$ being the stretch ratios along $\underline{N}_1$ and $\underline{N}_2$ directions respectively. Finally, the hyperelastic potential has been refined from Charmetant's initial proposal [16]. The revised expression Eq. (15) is designed to model an asymmetric (but continuously differentiable) response, thereby distinguishing between compaction and expansion. Indeed, the compaction is much stiffer than expansion, as all fibers are in contact in the former case, whereas no cohesive interaction (except for a slight effect of the binder) occurs in the latter case. The hyperelastic potential for compaction then writes

$$W_{\text{comp}} = \begin{cases} K_{\text{comp}} |I_{\text{comp}}|^p & \text{if } I_{\text{comp}} \le 0 \\ K_{\text{expan}} I_{\text{comp}}^p & \text{if } 0 < I_{\text{comp}} \le I_{\text{thold}} \\ \frac{1}{2} K_{\text{expan}} I_{\text{comp}}^2 & \text{if } I_{\text{thold}} < I_{\text{comp}} \end{cases} \qquad (15)$$

where $K_{\text{comp}}$ and $K_{\text{expan}}$ respectively designate the compaction and expansion stiffnesses, and $p$ is a real parameter controlling the non-linearity of the compaction, usually modeled with power laws. The $C^1$ differentiability in conjunction with convexity condition implies that $p > 2$. The expansion part is assumed to be linear and $K_{\text{expan}}$ depends on the characteristics of the yarn, such as the twist and the sizing of the yarn. The threshold value $I_{\text{thold}} = \frac{K_{\text{expan}} (p-2)^{-1}}{p K_{\text{comp}}}$ is defined as a transition value to ensure the hyperelastic potential is $C^1$. The transition between compaction and expansion occurs when the compaction invariant changes sign (i.e. at $I_{\text{comp}} = 0$).

### 2.4.2. Distortion

Distortion mode, represented in Fig. 6(b), is associated to a variation of the shape of the yarn's cross-section (at constant area). As a counterpart of the compaction mode, it can be expressed by the deviatoric part of the strain tensor reduced to its two-dimensional form in the ($\underline{N}_1$, $\underline{N}_2$) transverse plane. The matrix expression and modal decomposition of the intrinsic deformation gradient Eq. (5) and (12) yield

$$\underline{f}_{\text{dist}} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \frac{f_{11}}{f_{22}} & 0 \\ 0 & 0 & \frac{f_{22}}{f_{11}} \end{bmatrix}_{(\underline{M}, \underline{N}_1, \underline{N}_2)} \qquad (16)$$

for the distortion part, with associated invariant $I_{\text{dist}}$ defined by

$$I_{\text{dist}} = \frac{1}{2} \ln \left( \frac{f_{11}}{f_{22}} \right) \qquad (17)$$

The mechanical behavior for this mode is represented by the following power law

$$W_{\text{dist}} = K_{\text{dist}} I_{\text{dist}}^q \qquad (18)$$
with $K_{dist}$ the distortion stiffness and $q$ a nonlinear parameter.

Contrary to distortion potential proposed in Charmetant's work, we decide to model this behavior with a power law ($q \ge 2$). Indeed, fibers rearrangement, friction and yarn's sizing failure are micro phenomena more likely be modeled by a power law than a classic elastic linear one. One can prove with a first order Taylor expansion that the second Piola-Kirchhoff stress tensor does converge to zero if the invariant does. And with condition Eq. (4), we can deduce that for any deformation $F \simeq f_{11} \ge f_{22} \Rightarrow I_{dist} \ge 0$. The competition between compaction and distortion modes will determine the variation in shape and size of the yarn's cross-section, and can be illustrated through elementary transverse compression tests, as further discussed in Section 3.1.

### 2.4.3. Longitudinal shear

Longitudinal shear corresponds to the fibers relative slippage along the mean fiber direction. Following the same reasoning as previously for the compaction and distortion modes, the formulation for longitudinal shear proposed by [16] is retrieved from Eq. (19)-(21) below:

$$f_{\text{shear}} = \begin{bmatrix} 1 & \frac{f_{M1}}{f_M} & \frac{f_{M2}}{f_M} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (19)$$

$$I_{\text{shear}} = \sqrt{\frac{f_{M1}^2 + f_{M2}^2}{f_M^2}} \quad (20)$$

$$W_{\text{shear}} = \frac{1}{2} K_{\text{shear}} I_{\text{shear}}^2 \quad (21)$$

A homogeneous longitudinal shear is represented in Fig. 6(c), where $f_{M1}$ and $f_{M2}$ components give the rotation of the hexahedral non-fibrous faces with respect to the transverse plane. To facilitate comprehension, the representation has been simplified by presenting only the rotation associated with $f_{M2}$ and assuming that the fibers are not elongated, given that they are quasi-inextensible. This explains the value of $f_M = 1$.

Competition between longitudinal shearing and fiber bending has a major impact in yarn's curvature. Indeed, kinematic response of a bent yarn is fully determined by its fibers bending stiffness, their interaction such as slippage or cohesion and their quasi-inextensibility. Therefore, if longitudinal shear stiffness is too high, fiber are not able to slide on each other and the yarn's curvature is directly modified. This case will also be discussed further in Section 3.2.

### 2.5. Integration of each mode

Subsequent to the modal decomposition of energy and the construction of the solid-beam element, a specific integration is proposed for each mode. Given that each mode does not necessitate the same interpolation and is not encoded in the same element, this study proposes a tailored modal integration process. Two versions of the element were developed: one with full integration, and a second with selective reduced integration.

The full integration version is shown in Fig. 7(a), and encompasses height integration points in the hexahedral element, two for the bending on each fictive fiber, and one for the tension. Since elongation is assumed to be constant on each fictive fiber section, one integration point per truss element is enough. The Bernoulli beams framework used in NEM implies that the curvature is linear and thus prescribes two integration points per fictive fiber for bending.

A second version of the element can be used with partial reduced integration (see Fig. 7(b)) for the sake of computational efficiency. When performing reduced integration in the hexahedral element, classical solution would be to use a unique integration point. Here, we propose to reduce partially the integration in the longitudinal direction, and keep a full integration in the transverse plane; in fact, considering the transverse isotropy assumption for the yarn's behavior,

compaction, distortion and longitudinal shear only need be evaluated in the transverse plane.

In case of a reduced integrated hexahedral element, spurious modes that need to be stabilized might occur [43]. In the case of reduced integration (one integration point) twelve hourglass modes exist. By keeping four integration points, only eight of them remain.

Usually, stabilization methods like assumed natural strain [44,45] or enhanced assumed strain [46] are applied to prevent the occurrence of hourglass modes. For this fibrous solid-beam element, the remaining eight zero-energy modes, shown in Appendix B (Fig. B.18), can be stabilized by the bending mode, as it has been shown in previous research for in-plane bending in macroscopic forming simulation [34,41]; this specific aspect will be further discussed later.

## 3. Analyzing modal contribution through numerical simulations

This section is devoted to a discussion of the relative importance of deformation modes in the kinematics of the solid-beam element. This analysis is conducted through the presentation of numerical simulations, which have been selected to highlight the potential of the fibrous solid-beam element and its associated hyperelastic constitutive law. All simulations were performed by using an explicit time integration scheme based on [47], which is implemented in an in-house finite element code. During each simulation, the kinetic energy was carefully monitored to ensure we remain within the quasi-static assumption. In Fig. 8, we consider the following three numerical tests: (a) a boundary-free transverse compression test, (b) a cantilever bending test, and (c) buckling test. In each simulation, the mean fiber direction $M$ is aligned with the global $X$ direction. The simulation are velocity-driven, with $V_x$ (resp. $V_z$ and $V_y$) denoting the prescribed constant velocity on the boundary in the $X$ (resp. $Z$ and $Y$) direction.

The choice of material parameters is discussed for each simulation, based on kinematic results only (stresses are not considered at this stage). For the sake of simplicity, only the relative influence of the modes stiffnesses will be analyzed in the present study. The influence of the nonlinear parameters $p$ and $q$ (see Eq. (15) and (18)) are easier to predict and of lesser interest.<fnref n="1" /> Therefore, they will be fixed to $p = 2.5$ and $q = 2.3$ in the rest of the paper. Next, several sets of material parameters representative of different classes of fibrous materials are studied. In what follows, we introduce the concept of parameter ratios, and we discuss the choice of these ratios in line with material behavior. This will lead us to categorize families of fibrous materials according to these ratios.

### 3.1. Incompressibility ratio

As introduced in Section 2.4.2, the ratio between compaction and distortion stiffnesses can have a great influence on the element kinematics. This is demonstrated here with an elementary transverse compression test. Boundary conditions are described in Fig. 8(a). Since fictive fibers are aligned with direction $M$ and the longitudinal shearing stiffness is not null $W_{shear}$, the test is based on plain strain assumption in $(Y, Z)$ plane.

The simulation is conducted with various stiffness ratios. First, $R_i$ ratio is introduced, referred to as the incompressibility ratio:

$$R_i = \frac{K_{\text{comp}}}{K_{\text{dist}}} \quad (22)$$

<footnote n="1">¹ For instance, in the case of boundary-free transverse compression, the compaction invariant $I_{comp}$ is most likely bounded ($-1 < I_{comp} < 0$), since a value of $I_{comp} < -1$ would indicate a 87% reduction of the cross-sectional area, which is hardly ever observed in fiber yarns. Thus, when $p$ increases, the compaction energy decreases, hence spreading is reduced because the material is increasingly compressible (see Eq. (15)).</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Solid-beam element versions build on hexahedral solid element with additional fictive fibers on edges 1–2, 3–4, 5–6, 7–8: (a) full integration and (b) partial reduced integration. Integration points are specific to each mode: (■) integration in the volume, (▲) integration for bending, (●) integration for fiber tension.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Boundary conditions of the three simulations performed to analyze modal contribution to yarn behavior: (a) boundary-free transverse compression test, (b) cantilever bending test, (c) buckling test.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Influence of compressibility ratio $R_i$ on transverse compaction spreading (left) and deformed configurations (right), with (top) $R_i = 1$, i.e. $K_{comp} = K_{dist}$ and (bottom) $R_i = 500$ i.e. $K_{comp} = 500 K_{dist}$.</figcaption></figure>

The results are displayed in Fig. 9 by focusing on the element spreading: half of the initial state of the element is reproduced in solid black line while the compacted state is represented in dashed colored line for each pair of ($K_{comp}$, $K_{dist}$) parameter values corresponding to ratios $R_i = 1$ and $R_i = 500$. As the deformed states are almost superimposed, we provide the final width of the yarn for each case in the left bar charts.

Fig. 9 shows that spreading is almost independent of the stiffness value for a given ratio $R_i$. Note that the stresses, of course, would vary with both stiffnesses.

The same simulation is conducted to see the influence of the compaction material parameter $K_{comp}$ on the yarn’s kinematic response, with a fixed distortion stiffness $K_{dist} = 10^{-5}$. Fig. 10 highlights the competition between compaction and distortion modes for yarn transverse compaction: when $R_i$ increases, the behavior converges to nearly incompressible material. This is illustrated in Fig. 10(a), where the maximum spread reaches the theoretical limit of an incompressible

material for $R_i = 10^5$; at the value of $R_i = 10^3$, the spread is already at 98% of this incompressibility limit. This gives a clear indication of the range in which one can set the $R_i$ ratio depending on the fibrous material under consideration.

## 3.2. Longitudinal shear ratio

Fibrous stacks are known to exhibit specific decoupling between bending and longitudinal shear stiffness. While the bending stiffness $K_{bend}$ drives the mean deflection of a beam, the longitudinal shear stiffness plays an important role in the occurrence of fiber-to-fiber slippage. A longitudinal shear ratio $R_t$ is therefore introduced to control this effect.

$$R_t = \frac{K_{shear}}{K_{bend}} \qquad (23)$$
<figure><img src="image_11.png" /><figcaption>Fig. 10. Elementary simulation of compaction spreading with various incompressibility ratios $R_i$: (a) Normalized spreading against $R_i$ (a value of 1 means no spreading); theoretical limit for incompressible behavior (100%) and 98% of this value are plotted in dashed and dotted lines respectively; (b) initial and final configurations for various $R_i$ ratios.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Effect of longitudinal shear ratio $R_i$ in a cantilever bending test. (a) High ratios lead to Bernoulli-like kinematics, whereas (b) low ratios (negligible shear stiffness) lead to shear between fibers. The other parameters are fixed $C_{10} = K_{comp} = K_{dist} = K_{expan} = K_{shear}$.</figcaption></figure>

The role of $R_i$ is illustrated by a cantilever bending test. A velocity is applied on two nodes at the free end, whereas the four nodes at the other end are clamped, as shown in Fig. 8(b). The simulation is 3D, but the results are plotted in the $(Z, X)$ plane; in fact for the given set of parameters, displacements in $Y$ direction were negligible. In order to discuss the role of $R_i$, the other material parameters are kept equal: $C_{10} = K_{comp} = K_{dist} = K_{expan} = K_{shear}$.

The influence of $R_i$ on the cross-section orientation is visible in Fig. 11. The two conditions $R_i \gg 1$ and $R_i \ll 1$ led to two extreme cases: on the one hand, the cross-sections remain orthogonal to the neutral axis, following Bernoulli-like kinematics (Fig. 11(a)); on the other hand, the cross-sections are free to rotate (Fig. 11(b)), which is more representative of a pile of independent fibers with interlayer slippage, or a stack of paper leaves [48].

Note that in the case where $R_i \gg 1$, in order to enable top fiber elongation and bottom fiber compaction, $C_{10}$ had to be of the same order of magnitude as the other stiffnesses. Another solution might be to allow for variations in thickness, which will be discussed later.

### 3.3. Inextensibility ratio

In the case of technical textiles, fibers are often considered inextensible. The parameters set used above does not guarantee fiber quasi-inextensibility. In fact, as we move towards the Bernoulli-like kinematics of Fig. 11(a), the fictive upper and lower fibers are subjected to tension and compression respectively. A third ratio $R_X$ is therefore

introduced to drive the inextensibility of the fibers when needed. It is defined as

$$R_X = \frac{C_{10}}{K_{comp}} \qquad (24)$$

In the study, a ratio of $R_X = 100$ was shown to be sufficient to ensure quasi-inextensibility of the fibers. In the context of explicit dynamic schemes, it is indeed recommended not to increase this parameter excessively in order to maintain a reasonable critical time increment.

Fig. 12 shows the solution of a cantilever test with $R_X = 100$. The solution includes local buckling near the clamped edge. This local buckling is attributed to the incompatibility between fiber quasi-inextensibility, local curvature and absence of relative sliding. This is further discussed after introducing the last parameter ratio in the next section.

### 3.4. Buckling ratio

Previous simulations have shown that fiber inextensibility coupled with a high shear stiffness leads to local fiber buckling. Deformation of the beam's cross-section can also be a key element to reproduce fibrous kinematics. While compaction could be hard in certain cases (like in paper stacks), the expansion of the cross-section, i.e. an opening of a fiber stack can happen, in particular in fiber bundles without binder.

In order to control this buckling effect, a fourth ratio is introduced, called the buckling ratio $R_b$. This ratio can be split into two values, $R_b^{-}$
<figure><img src="image_13.png" /><figcaption>Fig. 12. Kinematic response on a cantilever bending test. Local fiber buckling occurs due to the incompatibility between fiber quasi-extensibility, local curvature and absence of relative sliding. Parameters are: $R_i = 1$, $R_t = 10^3$, $R_X = 100$.</figcaption></figure>

for compaction and $R_b^+$ for expansion, defined as

$$R_b^- = \frac{K_{\text{comp}}}{K_{\text{bend}}}, \qquad (25a)$$

$$R_b^+ = \frac{K_{\text{expan}}}{K_{\text{bend}}}. \qquad (25b)$$

The effect of $R_b$ is shown on a buckling test (see Figs. 8(c) and 13). Compaction stiffness is kept high, but expansion stiffness is modified through $R_b$. The lower the ratio, the lower is the cohesion between fibers; this allows buckling to appear at a structural scale, reproducing de-cohesion of the fiber bundle. This test may seem a little far from the use of fiber yarns in textile applications. However, we show next that this phenomenon has a determining impact on roving deformation, by returning to the cantilever bending test.

Fig. 14 shows again the results of the cantilever bending test with a variety of parameter sets. The cantilever test of Fig. 11 was performed again with other incompressibility ratio $R_i$ and longitudinal shear ratio $R_t$, and considering quasi-inextensible fibers ($R_X = 10^2$) and low fiber cohesion ($R_b^+ = 10^{-1}$, $R_b^- = 10^{-4}$). The results focus on the compaction and distortion invariants, i.e. the shape of the cross-section. Fig. 14(a) shows that section changes are inevitable, by allowing compressibility ($R_i = 1$) or allowing free sliding between incompressible fibers ($R_i = 10^{-4}$ and $R_i = 10^3$). Fig. 14(b) shows $I_{\text{dist}}$ values, focusing on the influence of $R_i$ on the cross-section variations. When $R_i$ is high (i.e. $K_{\text{comp}} \gg K_{\text{dist}}$), the cross-section cannot be compressed. Thus, the inextensibility of fibers lead to isochoric cross-section change, resulting from fiber local buckling.

These simulations show what happens in an inextensible fiber bundles: bending of the bundle (structural effect) will cause local buckling of fibers and section changes along the fiber yarns, leading to a heterogeneous fiber density field (as shown by $I_{\text{comp}}$ distribution on Fig. 14(a)).

### 3.5. Torsion and buckling of a fiber bundle

In this section, we now consider a common solicitation of fiber bundles in textiles which would involve a coupling of the different deformation modes: torsion and buckling. As these examples deal with auto contact, line-line contact elements are attached to the edges and node-surface contact elements are attached to the faces of the solid element.

#### 3.5.1. Torsion-buckling

The test shown on Fig. 15 is a two-step simulation starting with a torsion (step 1), followed by a buckling load (step 2). Boundary conditions are represented in Fig. 15(a). Material parameters chosen such that $R_i = 1$, $R_X = 10^2$, $R_t = 10^{-2}$, $R_b = 10^{-2}$.

During the torsion step (Fig. 15(b)), compaction and distortion stiffness are low enough to allow fictive fibers to be in contact with one another in the center part. During the buckling step (Fig. 15(c)),

the expansion and distortion stiffness are low enough to obtain de-cohesion in the center. Invariant $I_{\text{comp}}$ is plotted on both figures to highlight changes in the cross-sectional area, with positive (respectively negative) values indicating expansion (respectively compaction).

#### 3.5.2. Buckling-torsion

A second test is made, inverting the two steps. Now the buckling step is performed first, then the torsion follows. This test is equivalent to the torsional post-buckling of elastic rods presented in [49].

On Fig. 16 we show two steps of this simulation: the state after buckling, and the final post-buckling twisting. One can observe the typical loop formation and auto-contact of the fiber bundle. Even with a coarse mesh, the solution is well reproduced by the proposed model.

### 3.6. Textile compaction

The last simulation example deals with the case of a sample of plain weave textile made of 8 weft yarns and 8 warp yarns, which is compressed between two rigid surfaces (see Fig. 17). In a first step, starting from perfectly flat yarns, we mimic the weaving process to obtain a textile with imperfect geometry where yarn's density and crimp are intentionally heterogeneous. Thus, this sample has no plane of symmetry. It should be noted that the mesh generation method used herein has enabled us to generate a textile with an initial state of stress. The next simulation step is the compaction of the textile between the two rigid surfaces.

During the compression step, every node is totally free, and yarns are free to slide on each other, with a yarn-yarn and yarn-rigid surface frictionless condition. The absence of friction is not an hypothesis that fits real applications, especially in the case of yarn-yarn friction. In fact, yarn-yarn friction has been demonstrated to be a multiscale phenomenon (micro-scale fibers, and mesoscale yarns), and to be highly anisotropic [50]. The mechanism is even more complex when considering twisting of yarns [51]. Considering that modeling properly this kind of contact is still a challenging task, most of the existing simulations simplify contact by considering it homogeneous and isotropic. In order to better focus our attention on the behavior law and stability of our element, we therefore chose to perform simulations without friction.

This simulation was run with an incompressibility ratio $R_i = 10^3$. The compaction invariant field is plotted on the deformed configuration in Fig. 17(b). Invariant $I_{\text{comp}}$ Eq. (14) was chosen because it quantifies the cross-sectional area ratio between the initial and current configurations, which could be linked to the variation of fiber volume fraction. As can be seen from Fig. 17(b), the lack of symmetry in the initial geometry resulted in a non-homogeneous compaction state in the textile. This simulation demonstrates, as a proof-of concept, the capability of the solid-beam element to capture local variations in fiber volume fraction resulting from the forming process.

## 4. Discussion

The solid-beam element proposed here is capable of modeling different types of fiber stacking. Depending on the value of property ratios, it can represent bundles with weak cohesion (as a strand of hairs) or with a bit more cohesion (like carbon yarns with binder). The cohesion can be anisotropic, splitting the transverse section cohesion (driven by the compaction mode), and the longitudinal shearing cohesion (driven by longitudinal shearing mode). The role of longitudinal shear stiffness was discussed previously in the literature, and experimental tests were proposed (see [18,52]), showing that the material parameter is rather low; although this deformation mode is rarely used intentionally in weaving or shaping processes, it can appear during the process and play an important role in the degradation of roving over time [51,52].

The de-cohesion in the transverse section is also a mode that should be considered in simulations. This is made possible with our element, and dedicated material parameters. Several simulations presented in
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Buckling test with various buckling ratios on a fiber bundle. (a) $R_b = 1$; (b) $R_b = 10^{-4}$. Color bar shows compaction invariant $I_{comp}$. Note that the yarn's cross-section is expanding as $I_{comp} > 0$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Cantilever bending test on a fiber bundle with various parameter sets. Constant parameters are $R_X = 10^2$, $R_b^+ = 10^{-1}$, $R_b^- = 10^{-4}$: (a) Deformed shape showing compaction invariant $I_{comp}$; (b) zoom for $R_t = 1$, showing distortion invariant $I_{dist}$.</figcaption></figure>

this work have shown that fiber de-cohesion in the transverse plane is related to fiber bending phenomenon and incompatibility with their inextensibility. From this point of view, the model provides a better understanding of the effects observed experimentally, for instance the changes in cross-section showed in [51]. Looking back to Fig. 14, we see that the cross-section shape modifications can be quite complex depending on boundary conditions and material parameters. Fiber distribution and density during forming process could be strongly affected by these phenomena.

Another strong point of the proposed model is that it is able to consider different torsional states for a bundle of fibers. This allows to model the effect of initial torsion on the mesoscopic behavior of the yarn alone or embedded in a reinforcement.

A set of parameter ratios were proposed, and were shown to drive the response of a fiber bundle. Table 1 gives the lower and higher values of ratios that have been used in this paper. The table also summarizes the main effect of a given ratio on the fiber bundle behavior. This table can be used in modeling context, for instance when using the proposed approach, but could also be used as guidelines when performing experimental identification.

Two versions of the element were proposed: one with full integration in the solid element, and one with partial reduced integration. Concerning the partial (longitudinal) reduced integration, adding fictive fiber was shown to stabilize zero energy modes; this principle

had already been evoked in adding in-plane bending in shells [41]. However, one should keep in mind that this beneficial side effect only solves the problem of zero energy modes. The equivalent hourglass modes can still be activated when energy balance allows it. In this configuration, bending stiffness will act as a penalty coefficient to the solid hourglass mode, which means that the bending stiffness cannot be selected irrespective of the other stiffnesses. Depending on the boundary conditions, one could need to stiffen the bending behavior to that purpose. In general, we thus recommend the complete integration to avoid this problem. For a simulation such as the bending cantilever test, the model based on partially reduced integrated element offered a 30% faster resolution. Note that the computations was performed without improving efficiency (e.g. with parallelization). Considering all these points, the authors recommend using a fully integrated version of the element. Using a partially reduced integrated version was however helpful to understand and detail the contribution of the bending mode to the mechanics of our element.

The new solid-beam element was used to model a portion of a textile (Fig. 17); in the context of this paper, this application was limited to a small sample. Apart from computation time, there are no obstacles to apply it to larger samples, or even to multi-layer forming. However, the model proposed here assumes the existence of yarns, and may exclude certain non-woven reinforcements, like uni-directional reinforcements, for which specific solid-shell approaches could be developed.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Sequential torsion-buckling simulation on a fiber bundle: (a) Boundary conditions for step 1 and 2; (b) deformed shape after step 1 and (c) after step 2. Color bar shows the compaction invariant.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Torsional post-buckling simulation. Initial configuration and mesh are represented in wireframe. Two steps of the simulation are superimposed: buckling step and final torsional post-buckling state. The color bar shows distortion invariant value for both states.</figcaption></figure>

<table><caption>**Table 1**

Synthesis on material properties ratios.</caption><thead><tr><th>Ratio</th><th>Role</th><th>Low</th><th>Effect</th><th>High</th><th>Effect</th></tr></thead><tbody><tr><td>R<sub>x</sub></td><td>fiber inextensibility</td><td>1</td><td>stretched fibers</td><td>10<sup>3</sup></td><td>inextensible</td></tr><tr><td>R<sub>i</sub></td><td>section incompressibility</td><td>1</td><td>compressible</td><td>10<sup>5</sup></td><td>incompressible</td></tr><tr><td>R<sub>t</sub></td><td>longitudinal shear</td><td>10<sup>-4</sup></td><td>free sliding</td><td>1</td><td>no sliding</td></tr><tr><td>R<sub>b</sub><sup>+</sup></td><td>expansion buckling</td><td>10<sup>-4</sup></td><td>de-cohesion possible</td><td>1</td><td>no de-cohesion</td></tr><tr><td>R<sub>b</sub><sup>-</sup></td><td>compaction buckling</td><td>10<sup>-3</sup></td><td>collapsing</td><td>1</td><td>encourage yarn buckling</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Unit cell of a plane weave textile simulations, made of 8 weft yarns and 8 warp yarns; (a) weaving process and (b) textile after compaction with incompressible yarns (top view in (X, Y) plane). Color shows the compaction invariant.</figcaption></figure>

## 5. Conclusion

A new solid-beam dedicated to fiber yarns was proposed. It consists in an 8-nodes hexahedral element assembled with four fictive fibers on the edges along the fiber direction. Material behavior is based on transversely isotropic hyperelastic law within the solid element, completed with fiber bending and tension energy on the fictive fiber edges. The category of fibrous element introduced in [53] is thus extended to include a solid-beam element. The approach using fictive fibers and the neighboring elements method to take into account the bending energy in the fibrous medium has proved to be a relevant alternative to second gradient models.

Working a new finite element, this study led us to improve the formulation of hyperelastic media with fiber directions, and to analyze the role of material parameters in fibrous materials response. To do so, material parameter ratios were introduced, and their influence on the fibrous behavior was discussed. This original approach provides a better understanding of what distinguishes fiber bundles of different natures (more or less bounded). It could serve as a guide for both simulation and experimental material characterization.

The element shows its ability to describe fibrous beam in a great variety of applications, including simulation of woven reinforcements. In future work, this approach can be used to perform forming simulation at the macroscopic level in order to obtain realistic maps of fiber density.

## CRediT authorship contribution statement

**Baptiste Lacroix:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Julien Colmars:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Conceptualization. **Auriane Platzer:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Methodology, Conceptualization. **Emmanuelle Vidal-Sallé:** Writing – review & editing, Supervision.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Baptiste Lacroix reports financial support was provided by French National Research Agency. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This research was funded in part by the French Agence Nationale de la Recherche, ANR, project ANR-21-CE06-0040-01 COQFIB. A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions.

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

## Appendix A. About local basis in transverse section

First, let us consider a local orthonormal basis ($M, Z_1, Z_2$), attached to the undeformed initial state, with $M$ the local fiber mean direction. Vectors ($Z_1, Z_2$) can have any orientation, thus we chose to build them by using global orthonormal referential $(x, y, z)$, see Eq. (A.1).

$$ \begin{aligned} \underline{Z}_2 &= \underline{M} \times (\underline{z} \times \underline{M}) \\ \underline{Z}_1 &= \underline{Z}_2 \times \underline{M} \end{aligned} \quad (A.1) $$

Note that the case where $M$ and $z$ are aligned is treated separately.

To get ($N_1, N_2$) we must first calculate the angle $\omega$ Eq. (A.2), which gives the rotation compared with ($Z_1, Z_2$). In this equation $C^{-1}$ is the inverse of right Cauchy-Green tensor.

$$ 2\omega = \operatorname{atan2}(-2\underline{Z}_2 \cdot \underline{C}^{-1} \cdot \underline{Z}_1, \underline{Z}_2 \cdot \underline{C}^{-1} \cdot \underline{Z}_2 - \underline{Z}_1 \cdot \underline{C}^{-1} \cdot \underline{Z}_1) \quad (A.2) $$

Finally, we can deduce ($N_1, N_2$) thanks to Eq. (A.3).

$$ \begin{aligned} \underline{N}_1 &= \cos \omega \underline{Z}_1 + \sin \omega \underline{Z}_2 \\ \underline{N}_2 &= -\sin \omega \underline{Z}_1 + \cos \omega \underline{Z}_2 \end{aligned} \quad (A.3) $$

Every component of the deformation gradient tensor expressed in ($M, N_1, N_2$), see Eq. (5), is then obtained by projecting the tensor on the corresponding dyadic product of the basis vectors.

## Appendix B. Zero energy modes in reduced integration

The element can be used in a full integration version, or a partially (along fiber direction) reduced integration version. In the latter case, eight zero energy modes can occur in the hexahedral element, as shown in Fig. B.18. However, in the context of the solid-beam element, the bending of fictive fibers will stabilize these modes.
<figure><img src="image_20.png" /><figcaption>Fig. B.18. Hourglass modes in case of partial reduced integration (along fiber direction). These modes can be stabilized by bending and tension modes adapted from [54].</figcaption></figure>

## Data availability

Data will be made available on request.

## References

[1] Sozer E, Simacek P, Advani S. Resin transfer molding (RTM) in polymer matrix composites. In: Manufacturing techniques for polymer matrix composites (pMCs). Elsevier; 2012, p. 245-309. http://dx.doi.org/10.1533/9780857096258.3.243, URL https://linkinghub.elsevier.com/retrieve/pii/B9780857090676500092.

[2] Arbter R, Beraud J, Binetruy C, Bizet L, Bréard J, Comas-Cardona S, Demaria C, Endruweit A, Ermanni P, Gommer F, Hasanovic S, Henrat P, Klunker F, Laine B, Lavanchy S, Lomov S, Long A, Michaud V, Morren G, Ruiz E, Sol H, Trochu F, Verleye B, Wietgrefe M, Wu W, Ziegmann G. Experimental determination of the permeability of textiles: A benchmark exercise. Compos Part A: Appl Sci Manuf 2011;42(9):1157-68. http://dx.doi.org/10.1016/j.compositesa.2011.04.021, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X11001357.

[3] Caglar B, Orgéas L, Rolland Du Roscoat S, Sozer EM, Michaud V. Permeability of textile fabrics with spherical inclusions. Compos Part A: Appl Sci Manuf 2017;99:1-14. http://dx.doi.org/10.1016/j.compositesa.2017.03.031, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X17301343.

[4] Cataldi M, Wielhorski Y, Moulin N, Parret-Fréaud A, Pucci MF, Liotier P-J. Accounting for mesoscale geometry and intra-yarn fiber volume fraction distribution on 3D angle-interlock fabric permeability. Int J Multiph Flow 2024;173:104721. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2024.104721, URL https://www.sciencedirect.com/science/article/pii/S030193222400003X.

[5] Wielhorski Y, Mendoza A, Rubino M, Roux S. Numerical modeling of 3D woven composite reinforcements: A review. Compos Part A: Appl Sci Manuf 2022;154:106729. http://dx.doi.org/10.1016/j.compositesa.2021.106729, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X21004437.

[6] Durville D, Baydoun I, Moustacas H, Périé G, Wielhorski Y. Determining the initial configuration and characterizing the mechanical properties of 3D angle-interlock fabrics using finite element simulation. Int J Solids Struct 2018;154:97-103. http://dx.doi.org/10.1016/j.ijsolstr.2017.06.026, Multi-scale Modelling of Fibrous and Textile Materials, URL https://www.sciencedirect.com/science/article/pii/S0020768317302925.

[7] Daelemans L, Tomme B, Caglar B, Michaud V, Van Stappen J, Cnudde V, Boone M, Van Paepegem W. Kinematic and mechanical response of dry woven fabrics in through-thickness compression: Virtual fiber modeling with mesh overlay technique and experimental validation. Compos Sci Technol 2021;207:108706. http://dx.doi.org/10.1016/j.compscitech.2021.108706, URL https://www.sciencedirect.com/science/article/pii/S0266353821000622.

[8] Simon J, Hamila N, Comas-Cardona S, Binetruy C. A yarn-scale woven fabric model including significant slippage formulated within the Arbitrary Lagrangian Eulerian framework. Comput Methods Appl Mech Engrg 2024;420:116741. http://dx.doi.org/10.1016/j.cma.2023.116741, URL https://linkinghub.elsevier.com/retrieve/pii/S0045782523008642.

[9] Iwata A, Inoue T, Naouar N, Boisse P, Lomov SV. Coupled meso-macro simulation of woven fabric local deformation during draping. Compos Part A: Appl Sci Manuf 2019;118:267-80. http://dx.doi.org/10.1016/j.compositesa.2019.01.004, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X19300041.

[10] Nguyen Q, Vidal-Sallé E, Boisse P, Park C, Saouab A, Bréard J, Hivet G. Meso-scopic scale analyses of textile composite reinforcement compaction. Compos Part B: Eng 2013;44(1):231-41. http://dx.doi.org/10.1016/j.compositesb.2012.05.028, URL https://linkinghub.elsevier.com/retrieve/pii/S135983681200354X.

[11] Wang J, Wang P, Hamila N, Boisse P. Mesoscopic analyses of the draping of 3D woven composite reinforcements based on macroscopic simulations. Compos Struct 2020;250:112602. http://dx.doi.org/10.1016/j.compstruct.2020.112602, URL https://linkinghub.elsevier.com/retrieve/pii/S0263822320311314.

[12] Hivet G, Boisse P. Consistent 3D geometrical model of fabric elementary cell. Application to a meshing preprocessor for 3D finite element analysis. Finite Elem Anal Des 2005;42(1):25-49. http://dx.doi.org/10.1016/j.finel.2005.05.001, URL https://linkinghub.elsevier.com/retrieve/pii/S0168874X05000600.

[13] Naouar N, Vidal-Sallé E, Schneider J, Maire E, Boisse P. Meso-scale FE analyses of textile composite reinforcement deformation based on X-ray computed tomography. Compos Struct 2014;116:165-76. http://dx.doi.org/10.1016/j.compstruct.2014.04.026, URL https://linkinghub.elsevier.com/retrieve/pii/S0263822314002189.

[14] Wijaya W, Ali M, Umer R, Khan K, Kelly P, Bickerton S. An automatic methodology to CT-scans of 2D woven textile fabrics to structured finite element and voxel meshes. Compos Part A: Appl Sci Manuf 2019;125:105561. http://dx.doi.org/10.1016/j.compositesa.2019.105561, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X19303100.

[15] Wijaya W, Kelly P, Bickerton S. A novel methodology to construct periodic multi-layer 2D woven unit cells with random nesting configurations directly from μ CT-scans. Compos Sci Technol 2020;193:108125. http://dx.doi.org/10.1016/j.compscitech.2020.108125, URL https://linkinghub.elsevier.com/retrieve/pii/S0266353819330453.

[16] Charmetant A, Vidal-Sallé E, Boisse P. Hyperelastic modelling for mesoscopic analyses of composite reinforcements. Compos Sci Technol 2011;71(14):1623-31. http://dx.doi.org/10.1016/j.compscitech.2011.07.004, URL https://linkinghub.elsevier.com/retrieve/pii/S0266353811002351.

[17] Criscione JC, Douglas AS, Hunter WC. Physically based strain invariant set for materials exhibiting transversely isotropic behavior. J Mech Phys Solids 2001;49(4):871-97. http://dx.doi.org/10.1016/S0022-5096(00)00047-8, URL https://linkinghub.elsevier.com/retrieve/pii/S002250960000478.

[18] Zheng R, Naouar N, Colmars J, Platzer A, Schäfer B, Morestin F, Kärger L, Boisse P. Mesoscopic finite element modeling of biaxial non-crimp fabric including representative stitch pattern. Compos Struct 2024;339:118126. http://dx.doi.org/10.1016/j.compstruct.2024.118126, URL https://linkinghub.elsevier.com/retrieve/pii/S026382232400254X.

[19] Gatouillat S, Bareggi A, Vidal-Sallé E, Boisse P. Meso modelling for composite preform shaping – Simulation of the loss of cohesion of the woven fibre network. Compos Part A: Appl Sci Manuf 2013;54:135-44. http://dx.doi.org/10.1016/j.compositesa.2013.07.010, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X13001942.

[20] Wang Y, Sun X. Digital-element simulation of textile processes. Compos Sci Technol 2001;61(2):311-9. http://dx.doi.org/10.1016/S0266-3538(00)00223-2, URL https://linkinghub.elsevier.com/retrieve/pii/S0266353800002232.

[21] De Luycker E, Hamila N. Modeling of hyperelastic bending of fibrous media using second-gradient isogeometric analysis: Weaving and braiding applications. Compos Part A: Appl Sci Manuf 2023;167:107415. http://dx.doi.org/10.1016/j.compositesa.2022.107415, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X22005966.
[22] 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, URL https://linkinghub.elsevier.com/retrieve/pii/S0020768316303936.

[23] Moustacas H, Durville D, Wielhorski Y. Modélisation et simulation par élément fini du comportement transverse de mèches de fibres de carbone. In: 13e Colloque national en calcul des structures, Université Paris-Saclay, May 2017, Giens, Var, France. 2017.

[24] Moustacas H, Durville D, Wielhorski Y. Enrichissement d'une cinématique poutre. Applications aux textiles en carbone.. In: 14ème Colloque National en Calcul de Structures (CSMA 2019), CSMA, LEM3, MSME, Université de Lorraine, Arts et Métiers, CNRS, May 2019, Hyères, France. 2019.

[25] Hu Y, Zhao Y, Liang H. Refined beam theory for geometrically nonlinear Pre-twisted structures. Aerospace 2022;9(7):360. http://dx.doi.org/10.3390/aerospace9070360, URL https://www.mdpi.com/2226-4310/9/7/360.

[26] Frischkorn J, Reese S. A solid-beam finite element and non-linear constitutive modelling. Comput Methods Appl Mech Engrg 2013;265:195-212. http://dx.doi.org/10.1016/j.cma.2013.06.009, URL https://linkinghub.elsevier.com/retrieve/pii/S0045782513001618.

[27] Wei G, Lardeur P, Druesne F. A new solid-beam approach based on first or higher-order beam theories for finite element analysis of thin to thick structures. Finite Elem Anal Des 2022;200:103655. http://dx.doi.org/10.1016/j.finel.2021.103655, URL https://linkinghub.elsevier.com/retrieve/pii/S0168874X21001360.

[28] Xiong H, Guzman Maldonado E, Hamila N, Boisse P. A prismatic solid-shell finite element based on a DKT approach with efficient calculation of through the thickness deformation. Finite Elem Anal Des 2018;151:18-33. http://dx.doi.org/10.1016/j.finel.2018.08.003, URL https://linkinghub.elsevier.com/retrieve/pii/S0168874X18302968.

[29] Betsch P, Stein E. An assumed strain approach avoiding artificial thickness straining for a non-linear 4-node shell element. Commun Numer Methods Eng 1995;11(11):899-909. http://dx.doi.org/10.1002/cnm.1640111104, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.1640111104, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.1640111104.

[30] Schwarze M, Reese S. A reduced integration solid-shell finite element based on the EAS and the ANS concept—Large deformation problems. Numer Meth Eng 2011;85(3):289-329. http://dx.doi.org/10.1002/nme.2966.

[31] Gao S, Liang B, Vidal-Salle E. Development of a new 3D beam element with section changes: The first step for large scale textile modelling. Finite Elem Anal Des 2015;104:80-8. http://dx.doi.org/10.1016/j.finel.2015.06.003, URL https://linkinghub.elsevier.com/retrieve/pii/S0168874X15000906.

[32] Ménard F, Cartraud P. Solid and 3D beam finite element models for the nonlinear elastic analysis of helical strands within a computational homogenization framework. Comput Struct 2021;257:106675. http://dx.doi.org/10.1016/j.compstruc.2021.106675, URL https://linkinghub.elsevier.com/retrieve/pii/S0045794921001978.

[33] Courteau-Godmaire H. A cosserat approach to modelling the deformation mechanisms of uncured composites during processing [Ph.D. thesis], Electronic theses and dissertations (eTDs) 2008+, University of British Columbia; 2023. http://dx.doi.org/10.14288/1.0436915, URL https://open.library.ubc.ca/collections/ubctheses/24/items/1.0436915.

[34] 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:111014. http://dx.doi.org/10.1016/j.ijsolstr.2021.03.001, URL https://linkinghub.elsevier.com/retrieve/pii/S0020768321000858.

[35] Liang B, Colmars J, Boisse P. A shell formulation for fibrous reinforcement forming simulations. Compos Part A: Appl Sci Manuf 2017;100:81-96. http://dx.doi.org/10.1016/j.compositesita.2017.04.024, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X17301744.

[36] Bai R, Colmars J, Naouar N, Boisse P. A specific 3D shell approach for textile composite reinforcements under large deformation. Compos Part A: Appl Sci Manuf 2020;139:106135. http://dx.doi.org/10.1016/j.compositesita.2020.106135, URL https://linkinghub.elsevier.com/retrieve/pii/S1359835X20303742.

[37] Badel P, Vidal-Sallé E, Maire E, Boisse P. Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale. Compos Sci Technol 2008;68(12):2433-40. http://dx.doi.org/10.1016/j.compscitech.2008.04.038, URL https://linkinghub.elsevier.com/retrieve/pii/S0266353808001486.

[38] Oñate E, Cervera M. Derivation of thin plate bending elements with one degree of freedom per node: a simple three node triangle. Eng Comput 1993;10(6):543-61. http://dx.doi.org/10.1108/eb023924, URL https://www.emerald.com/insight/content/doi/10.1108/eb023924/full/html.

[39] Brunet M, Sabourin F. A simplified triangular shell element with a necking criterion for 3-D sheet-forming analysis. J Mater Process Technol 1995;50(1-4):238-51. http://dx.doi.org/10.1016/0924-0136(94)01384-D, URL https://linkinghub.elsevier.com/retrieve/pii/092401369401384D.

[40] Brunet M, Sabourin F. Analysis of a rotation-free 4-node shell element. Internat J Numer Methods Engrg 2006;66(9):1483-510. http://dx.doi.org/10.1002/nme.1608, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.1608.

[41] Bai R, Guzman-Maldonado E, Zheng R, Colmars J. Influence of in-plane bending behaviour on textile composite reinforcement forming. Int J Mech Sci 2024;273:109206. http://dx.doi.org/10.1016/j.ijmecsci.2024.109206, URL https://linkinghub.elsevier.com/retrieve/pii/S0020740324002480.

[42] Bonet J, Wood RD. Nonlinear continuum mechanics for finite element analysis. Cambridge ; New York, NY, USA: Cambridge University Press; 1997.

[43] Belytschko T, Ong JS-J, Liu WK, Kennedy JM. Hourglass control in linear and nonlinear problems. Comput Methods Appl Mech Engrg 1984;43(3):251-76. http://dx.doi.org/10.1016/0045-7825(84)90067-7, URL https://linkinghub.elsevier.com/retrieve/pii/0045782584900677.

[44] Belytschko T, Bindeman LP. Assumed strain stabilization of the eight node hexahedral element. Comput Methods Appl Mech Engrg 1993;105(2):225-60. http://dx.doi.org/10.1016/0045-7825(93)90124-G, URL https://linkinghub.elsevier.com/retrieve/pii/004578259390124G.

[45] Abed-Meraim F, Combescure A. Stabilisation et projection de type " Assumed Strain " de l'élément de coque massif " SHBPS " pour éliminer les problèmes de verrouillage. In: 7e Colloque national en calcul des structures, May 2005, Giens, Var, France. 2005.

[46] Cardoso RPR, Yoon JW, Mahardika M, Choudhry S, Alves De Sousa RJ, Fontes Valente RA. Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements. Numer Meth Eng 2008;75(2):156-87. http://dx.doi.org/10.1002/nme.2250, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.2250.

[47] Belytschko T, Liu WK, Moran B, Elkhodary K. Finite elements for nonlinear continua and structures. John Wiley & Sons; 2014, URL https://scholar.google.com/scholar_lookup?title=Nonlinear%20finite%20elements%20for%20continua%20and%20structures&author=T.%20Belytschko&publication_year=2013.

[48] Huang J, Boisse P, Hamila N, Zhu Y. Simulation of wrinkling during bending of composite reinforcement laminates. Materials 2020;13(10):2374. http://dx.doi.org/10.3390/ma13102374, URL https://www.mdpi.com/1996-1944/13/10/2374.

[49] Thompson JMT, Champneys AR. From helix to localized writhing in the torsional post-buckling of elastic rods. Proc R Soc Lond Ser A Math Phys Eng Sci 1996;452(1944):117-38. http://dx.doi.org/10.1098/rspa.1996.0007, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1996.0007, URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1996.0007.

[50] Tourlonias M, Bueno M-A, Fassi G, Aktas I, Wielhorski Y. Influence of friction angle between carbon single fibres and tows: Experimental analysis and analytical model. Compos Part A: Appl Sci Manuf 2019;124:105478. http://dx.doi.org/10.1016/j.compositesita.2019.105478, URL https://www.sciencedirect.com/science/article/pii/S1359835X19302271.

[51] Walther J, Tourlonias M, Decrette M, Bueno M-A. Influence of multifilament yarn twist on yarn-to-yarn friction behaviour: Application to carbon fibre weaving. Compos Part A: Appl Sci Manuf 2023;174:107737. http://dx.doi.org/10.1016/j.compositesita.2023.107737, URL https://www.sciencedirect.com/science/article/pii/S1359835X23003135.

[52] Wendling-Hivet A, Ferré MR, Allaoui S, Nunez R, Loison S, Hivet G. Study of the cohesion of carbon fiber yarns: in-plane shear behavior. Int J Mater Form 2017-10-01;10(5):671-83. http://dx.doi.org/10.1007/s12289-016-1310-y, URL DOI: 10.1007/s12289-016-1310-y.

[53] Bai R, Colmars J, Chen B, Naouar N, Boisse P. The fibrous shell approach for the simulation of composite draping with a relevant orientation of the normals. Compos Struct 2022;285:115202. http://dx.doi.org/10.1016/j.compstruct.2022.115202, URL https://linkinghub.elsevier.com/retrieve/pii/S026382232200174.

[54] Abed-Meraim F, Combescure A. An improved assumed strain solid-shell element formulation with physical stabilization for geometric non-linear applications and elastic-plastic stability analysis. Numer Meth Eng 2009;80(13):1640-86. http://dx.doi.org/10.1002/nme.2676, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.2676.