Effects of Fiber Shape on Mechanical Properties of Fiber Assemblies

The effects of fiber shape on the mechanical responses of fiber assemblies under compression, tension, and shear deformations are numerically investigated using the discrete element method (DEM). Simulations of the compression of ring-shaped fibers are consistent with experimental results, verifying the discrete element method code. In the compressive tests of S-shaped fibers, pressure exhibits a nonmonotonic dependence on fiber curvature; while in the tensile tests, yield tensile stress generally decreases with increasing fiber curvature. In the shear tests, yield shear stress decreases with increasing fiber curvature for the S-shaped fibers, and the smallest yield shear stresses and the smallest coordination numbers are obtained for U-shaped and Z-shaped fibers. It is interesting to observe that for the assemblies of various fiber shapes, yield shear stress increases with increasing maximum Feret diameter of the fibers, which characterizes the largest dimension of a fiber between two parallel tangential lines. These novel observations of the effects of fiber shape provide some guidelines for material designs with the fibers.


Introduction
Fibrous granular materials are ubiquitous in industries of biomass fuels, agricultural products, polymers, textile products, fiber-reinforced composites, filter, and papers. Mechanical responses of fiber assemblies subject to external loads are critical for processing the materials and final outcomes of the products. Thus, previous studies were conducted to understand compression, tension, and shear flows of flexible fibers. Based on the extensive experimental results of fiber compressions, van Wyk [1] and Toll [2] proposed that applied pressure, P, and solid volume fraction, φ, of a fiber bed followed a power law relationship, and the exponent in the law, which is between 3 and 15.5, depended on the random or aligned microstructure of the fiber packings. Such power law relations of P and φ were also confirmed in the numerical simulations [3][4][5]. The simulation results showed that the increase in fiber-fiber friction coefficient caused the rapid rising of the pressure at a smaller solid volume fraction [4]. Coordination number, defined as an average number of contacting neighbors per fiber, and energy per fiber started to rise at a smaller solid volume fraction for the longer fibers than for the short ones [3]. Fiber stiffness had an impact on the compression behaviors [6][7][8][9]. The pressure curves as a function of sample strain became steeper as the fiber stiffness increased [6]. The addition of a small fraction of stiffer fibers to a bed of very flexible fibers could significantly increase the bulk stiffness of fiber bed by orders of magnitude [7]. For random networks of non-identical fibers, the overall network stiffness decreased as the variability of the fiber stiffness increased at constant mean fiber stiffness [8,9]. Smaller compressive loads and larger non-reversible deformation were obtained for plastic fibers rather than for elastic fibers, and it was found that fiber bending plasticity contributed to the smaller loads at early stage of compression and fiber-fiber contact plasticity caused smaller loads at late stage [10]. An increase in fiber-fiber contact stiffness resulted in a larger bulk stiffness of the fiber assembly [11], and the elasto-plastic contact force and adhesion synergistically enhanced the non-reversible deformation of the fiber bed after the compression [12]. In the numerical modeling of agricultural and biomass materials (e.g., crop stems), non-elastic, hysteretic contact forces, and fiber bending deformations should be considered for a more accurate prediction [13,14].
Some special fiber configurations were also considered in the previous studies. As reported by Masse and Poquillon [15], cross-link connections at the contact positions significantly increased the compression stresses and bulk stiffness at the early stage of compression, due to the strong constraints of cross-links on the relative motion between the connected fibers; nevertheless, larger compression stress and bulk stiffness were observed without cross-links at the late stage, which may be relevant to the stronger fiber orientational alignment without cross-links. Bergström et al. [16] and Borodulina et al. [17] observed that the tensile stiffness and tensile strength of a cross-linked fiber network depended on the fiber length distribution and material strength distribution. Persson and Isaksson [18] investigated dynamic fracture of cross-linked fiber networks under tensile loads, and they found that the continuum model was unavailable to describe the dynamic crack growth of the sparse fiber networks due to the very diffuse fracture zone. Li et al. [19] found that the addition of fibers to a spherical particle bed could increase shear strength of the granular material, enhancing the ability to bear larger loads, and a further improvement in the load bearing was achieved by arranging the fibers in a cage-like structure compared to in a random structure. Rodney et al. [20] performed both tension and compression tests on a packing of a self-entangled single long coiled wire, which showed unusual, reversible dilatancy behavior with the Poisson's ratio below 0 in tension and above 0.5 in compression.
The early work of Discrete Element Method (DEM) studies of flexible fibers included Yamamoto and Matsuoka [21] and Park and Kang [22]. The DEM method was further developed to investigate the compression of fiber assemblies [3][4][5][12][13][14], and microscopic fiber-scale information (fiber deformation, contact force, and local solid volume fraction) could be observed and analyzed. The DEM was also used to investigate shear flows of flexible fibers [23,24]. In dense flows at a constant solid volume fraction [23], steady shear stress increased with increasing friction coefficient, fiber stiffness, and fiber aspect ratio, as larger fiber-fiber contact forces were induced. In dense flows under a constant normal stress [24], the shear stress initially increased with the friction coefficient and then converged to an upper limit when the friction coefficient reached one, while the shear stress was insensitive to the fiber stiffness and fiber-fiber contact stiffness in such normal-stress-controlled flows.
In the previous studies, the fibers normally possessed the original shape of a straight line. The original shape is referred to as the equilibrium shape of a fiber under no external forces. How the shape of fibers affects mechanical behaviors of fiber assemblies is not yet well understood. However, a good understanding of the effects of fiber shape is crucial for material designs of the fibers, in order to achieve specific properties and functions. In the present work, the DEM method is used to simulate the compression, tension, and shear flows of the elastic fibers of various original shapes, including ring-shape, S-shape, U-shape, and Z-shape. The effects of the original shape of fibers on the bulk mechanical responses and microstructural properties of the fiber assemblies are analyzed based on the simulations results. The present article is organized as follows: The theoretical aspects and governing equations of the DEM-based flexible fiber model are presented in Section 2. The results of compressive tests, tensile tests, and shear tests are discussed in Section 3, Section 4, and Section 5, respectively. At last, the concluding remarks are drawn in Section 6.

Flexible Fiber Model
The flexible fiber model developed in previous work [6,10,25] has been used in the present study. As illustrated in Figure 1a, a fiber is discretized into several spherical nodes which are connected by cylinders of the same diameter to the sphere diameter. Virtual bonds are also used to connect the nodes, and they can deform as the fiber deforms. The deformation of bonds leads to the bond forces and moments, which act on the sphere nodes to resist the fiber deformation. The translational and rotational motion of a node sphere is governed by Newton's second law of motion: and in which v i and ω i are the translational and angular velocity vectors (m/s), respectively, of node sphere i with mass m i (kg) and moment of inertia J i (kg·m 2 ). The translational movement of the node sphere is driven by the normal contact force F c ni , tangential contact force F c ti , normal bond force F b ni , tangential bond force F b ti , contact damping forces F cd ni and F cd ti , bond damping forces F bd ni and F bd ti , and gravitational force m i g. Rotational movement is induced by the moments M c i , M b i , M cd i , and M bd i due to the contact forces, bond forces/moments, contact damping forces, and bond damping forces/moments, respectively. In the present fiber model, the basic element is a sphero-cylinder, as shown in Figure 1a. Thus, the contact detection between two fibers is based on the contact between two sphero-cylinders. As demonstrated in the previous work [10], three different contact types exist for a contact between two sphero-cylinders: hemisphere-hemisphere contact, hemisphere-cylinder contact, and cylinder-cylinder contact. The specific, geometry-dependent models of the normal and tangential contact forces, F c ni and F c ti , which were proposed in the previous work [10] to calculate the contact forces for each contact type, are utilized in the present simulations. The forces F on the right-hand side of Equation (1) have the unit of N, and the moments M on the right-hand side of Equation (2) have the unit of N·m. The bond forces and bond moments are functions of bond deformation, which is described by the relative displacements between two connected node spheres. Hence, the normal and tangential bond forces F b n and F b t are expressed as linear functions of normal and tangential displacements ∆ b n and ∆ b t , respectively, and The bond twisting moment M b twist and bond bending moment M b bend are computed incrementally based on the relative twisting angular velocity . θ twist and the relative bending angular velocity . θ bend between two connected node spheres and In where ν b is the Poisson's ratio of the bond) are the elastic modulus and shear modulus (Pa), respectively, of the bond material; A and l b are the cross-sectional area (m 2 ) and length (m), respectively, of the bond; I = πr 4 /4 is the area moment of inertia (m 4 ); I p = πr 4 /2 is the polar area moment of inertia (m 4 ); r is the radius of the fiber (m); and dt is the time step (s). An illustration of bond forces and moments acting on a node sphere is shown in Figure 1b.
Kinetic energy can be dissipated through deformation and vibration of the flexible fibers. This type of kinetic energy loss is implemented through bond damping forces and moments: and where K b n , K b t , K b twist , and K b bend represent the normal, shear, twisting, and bending stiffnesses (N/m), respectively, of the bond as defined in Equations (3) θ twist , and . θ bend , represent the relative normal velocity, tangential velocity, twisting angular velocity, and bending angular velocity, respectively, between two bonded node spheres of mass, m i and moment of inertia, J i . The kinetic energy dissipation rate due to the deformation and vibration of the flexible fibers is determined by the bond damping coefficient, β b . The larger β b , the faster the energy is dissipated.
The normal component F c n and tangential component F c t of the contact force F c exerted on a sphero-cylinder element are linearly distributed to the two node spheres of the element, as shown in Figure 1c. The normal and tangential components of the contact force on node sphere 1 can be expressed as and in which λ 1 and λ 2 are the distances between the contact point and the tangent points on the node spheres 1 and 2, respectively. Similarly, the force components on node sphere 2 have the expressions and The above algorithms were implemented in an in-house DEM code for the modeling of flexible fibers. To verify the code for the originally straight fibers, the uniaxial compression of a packing of cut rubber cords in a cylindrical container was simulated in the previous work [10], and the simulation and experimental results of loading/unloading curves are in good agreement. In the following section, the same DEM code will be validated for the modeling of originally crooked fibers.

Compression Tests
In this section, the curved fiber model is verified and calibrated against experimental results of compression tests of fiber rings. Thereafter, effect of fiber curvature on bulk compression behaviors is explored based on the simulations of S-shaped fibers.

Fiber Rings
To verify the DEM model of the fibers with curvature, experiments and simulations are performed on uniaxial compression of rubber rings. As shown in Figure 2a, a packed bed of 300 rubber rings is formed in a cylindrical container by releasing the rubber rings from the opening of the container one by one. The rubber rings have an outer diameter of circle of 23 mm and a diameter of line of 2.4 mm, a density of 1340 kg/m, and a Poisson's ratio of 0.5. The cylindrical container made of PMMA has an inner diameter of 80 mm and a height of 200 mm. After the initial packing, a steel punch, the top of which is connected to a ZwickRoell ® materials testing machine, moves down at a constant speed of 0.1 m/s into the container to compress the rubber ring bed, as shown in Figure 2b. The load exerted on the punch and the punch displacement are measured by the testing machine during the compression process. The same uniaxial compression of flexible rubber rings is simulated using the present DEM method, as shown in Figure 3, and the parameters used in the simulation are shown in Table 1. By performing a series of simulations with various fiber bending moduli E b , it is observed that the DEM simulation results are consistent with the experimental ones when the modulus E b is equal to 7.5 × 10 5 Pa. In the load-controlled compression, as shown in Figure 4a, the height of the granular bed decreases as the loading force increases, while in the displacement-controlled compression, as shown in Figure 4b, the pressure increases exponentially with the solid volume fraction, following a power law relationship like the originally straight fibers [1,2].

S-Shaped Fibers
As shown in Figure 5, S-shaped, U-shaped, and Z-shaped fibers are used in the present simulations to explore the effect of fiber shape on the mechanical characteristics of the fiber assemblies. In the simulations of the uniaxial compression of S-shaped fibers, as shown in Figure 6, the periodic boundary conditions are specified in the two horizontal directions (i.e., x and z directions), and the fibers are deposited under gravity on a flat wall. When the static state is achieved with negligible fiber velocities, the fiber bed is compressed by an upper wall, which moves downwards at a constant velocity of 0.1 m/s. It is noted that the compression results are insensitive to the loading speed in the range 0.01-0.1 m/s according to the present simulation results. The DEM simulation parameters are presented in Table 2. In the simulations, 300 S-shaped fibers of the linear length of 62.83 mm are generated in a rectangular domain of dimensions 60 × 300 × 60 mm 3 . The other properties of the S-shaped fibers are the same to those of the above rubber rings (see Table 1). The curvature of an S-shaped κ is defined as the reciprocal of the radius (i.e., 1/R), as shown in Figure 5.  The pressure P-solid volume fraction φ curves for the first three load-unload cycles of the straight fibers with κ = 0 m −1 and S-shaped fibers with κ = 100 m −1 are shown in Figure 7. The load-unload cycle 2 shifts to the right-hand side after the cycle 1 due to the fiber rearrangement and the consolidation of the fiber bed after the first load-unload cycle. The P − φ loop of the load-unload cycle 3 is very close to that of the cycle 2. At a given solid volume fraction φ, a larger pressure P is observed for the S-shaped fibers compared to the straight fibers, and a sharper increase in P with φ is obtained for the S-shaped fibers.

Parameters Values
Dimensions of domain, l x × l y × l z (mm 3 ) 60 × 300 × 60 Fiber shape S-shape, U-shape, Z-shape Linear length of a fiber, The P − φ curves in the first loading path and second loading path for the fibers with various curvatures κ are shown in Figure 8a,b, respectively. The straight fibers have the smallest pressures, and the loading curves of the S-shaped fibers exhibit a non-monotonic dependence on the fiber curvature. As shown in Figure 8c, the peak values of the pressure P at φ = 0.4 are obtained for the S-shaped fibers with κ = 44 m −1 and κ = 100 m −1 .
In the compression processes, coordination number, defined as the average number of contacting neighbors per fiber, increases with the solid volume fraction φ and pressure P, as shown in Figure 9a,b. At a specified φ, the coordination number in the load path is greater than that in the unload path (Figure 9a), due to the rearrangement of fibers. while at a specified P, the coordination number in the unload path is slightly greater than that in the load path (Figure 9b). Larger coordination numbers are obtained for the straight fibers with κ = 0 m −1 compared to the crooked fibers with κ = 100 m −1 . However, larger than average fiber-fiber contact forces are obtained for the crooked fibers with κ = 100 m −1 (Figure 10), causing larger pressures in the compression of the fibers with κ = 100 m −1 . Thus, the straight fibers have more but weaker contacts and hence lower pressures in the compression than the crooked fibers.  To examine the orientation of the fibers, an inclination angle of each fiber is measured as the angle of the line, which connects the two end points of a fiber, from the horizontal plane. The probability density distributions of the fiber inclination angles for the packings before and after the first load-unload cycle are shown in Figure 11a,b, respectively. Before loading, the fibers more likely have the inclination angles between 20 • and 60 • , while after the first load-unload cycle, the fibers are more likely oriented at the inclination angles between 10 • and 20 • . The curvatures κ of the fibers have limited impact on the distributions of the inclination angles. As shown in Figure 11c, the average inclination angles are significantly reduced after the first load-unload cycle, indicating that the fibers tend to align horizontally, which gives smaller potential energies of gravity. After the first load-unload cycle of compression, the fibers with the largest curvature considered (κ = 100 m −1 ) exhibit the largest average inclination angle.

Tensile Tests
The assembly of fibers is capable of bearing tensile loads due to strong interlocking and entanglement of the fibrous materials. Therefore, tensions of a bed of 500 fibers subject to stretching loads are also simulated. The DEM simulation parameters are the same to those used in the compression simulations, which are listed in Table 2. As shown in Figure 12, some fibers in the lower region (in gray color) are frozen and remain stationary. Some fibers in the upper region (in gray color), which are also frozen, move upwards in the y direction at a constant velocity of 0.1 m/s as a solid body to stretch the fibers in the middle region (in yellow color). In the tensile test simulations, the periodic boundary conditions are assigned in the horizontal x and z directions. Tensile stress σ yy is plotted against tensile strain ε yy in Figure 13a for the fibers with various curvatures κ. In general, the tensile stress initially increases and then decreases until the fiber assembly is torn apart. The peak value of the stress in the σ yy -ε yy curve is defined as the yield tensile stress σ 0 yy . After the peak σ 0 yy , the tensile stress σ yy decreases slowly with the strain ε yy for the fibers with κ ≤ 44 m −1 , and rapid decreases in σ yy are observed for the fibers with κ ≥ 58 m −1 . In addition, the yield tensile stress σ 0 yy generally decreases with increasing fiber curvature κ due to reduction in the interlocking of the fibers. The yield tensile stress σ 0 yy is dependent on the length L s /d and cross-sectional area A c of the cubic sample of fiber assembly. As shown in Figure 14, σ 0 yy increases with decreasing sample length L s /d and decreasing cross-sectional area A c : the smaller, the stronger. This size effect is due to the fact that more structural defects of weak fiber-fiber connections, which reduce the bulk strength, exist in a larger sample. It is also observed that the assemblies of straight fibers can hardly bear tensile stresses when the normalized sample length is larger than 15.

Shear Tests
Shear deformations of an assembly of 500 fibers with various fiber shapes are also simulated. The DEM simulation parameters are the same to those used in the compression and tensile simulations, which are listed in Table 2. As shown in Figure 15, some frozen fibers in the lower region (in gray color) remain stationary, and some frozen fibers in the upper region (in gray color) move horizontally in the x direction at a constant velocity of 0.1 m/s as a solid body to shear the fibers in the middle region (in yellow color). Periodic boundary conditions are used in the horizontal x and z directions. The evolution of shear stress σ xy with shear strain ε xy is plotted in Figure 16a for the straight fibers (κ = 0 m −1 ) with various sample lengths. The shear stress increases with the shear strain at the early stage and then fluctuates around a constant level, which is referred to as the steady state. At a given shear strain, the shear stress generally increases as the normalized sample length L s /d is reduced. In the shear process, the coordination number also increases with the shear strain at the early stage and then achieves a constant level, as shown in Figure 16b. Nevertheless, the coordination number, unlike the shear stress, shows a nonmonotonic and irregular dependence on the sample length.
Yield shear stress σ 0 xy , which is defined as the average shear stress at the steady state, is plotted in Figure 17a as a function of the normalized sample length L s /d for the straight fibers, S-shaped fibers with κ = 2, 44, and 100 m −1 , U-shaped fibers with θ = 0 • and 90 • , and Z-shaped fibers. Similar to the yield tensile stress (see Figure 14), the yield shear stress σ 0 xy decreases as the sample length L s /d increases. As shown in Figure 17a, for a given sample length, σ 0 xy decreases as the fiber curvature κ increases. The S-shaped fibers with a larger curvature κ = 100 m −1 , U-shaped fibers, and Z-shaped fibers have the smallest yield shear stresses. The average coordination numbers are insensitive to the normalized sample length, as shown in Figure 17b. In general, the coordination number decreases with increasing fiber curvature κ, and smaller coordination numbers are obtained for the U-shaped and Z-shaped fibers.  A Feret diameter, D F , is the distance between a pair of parallel lines which are tangential to the outline of a fiber particle, as shown in Figure 18a. By changing the direction of the parallel lines, different Feret diameters can be obtained for a specified fiber. Thus, a maximum Feret diameter, D max F , exists for a fiber, and it can vary as the fiber shape changes. In Figure 18b, the yield shear stress σ 0 xy is plotted against the normalized maximum Feret diameter, D max F /l f , in which the linear length of a fiber l f is equal to 62.83 mm for all the fibers in the shear deformation simulations. It can be seen that for the fibers of different shapes at various sample lengths, the yields shear stress tends to increase with increasing normalized maximum Feret diameter, D max F /l f . As a result, the shear properties of a fiber assembly are determined by the normalized maximum Feret diameter.

Conclusions
The numerical models of ring-shaped, S-shaped, and Z-shaped fibers are generated in the discrete element method (DEM) simulations. The mechanical responses of these fiber assemblies subject to compressive, tensile, and shearing loads are numerically investigated. The loading curves of bed height versus loading force (force loading) and pressure versus solid volume fraction (displacement loading) obtained from the compression simulations of rubber rings are consistent with the experimental results, verifying the developed DEM models for the crooked fibers.
In the compression tests, at a given solid volume fraction, higher pressures are obtained for the S-shaped fibers compared to the straight fibers (>5 kPa versus 4 kPa at solid volume fraction of φ = 0.4) due to the larger fiber-fiber contact forces. The effect of fiber curvature on the pressure is non-monotonic. The fiber assemblies can bear tensile forces due to the interlocking and entanglement of the fibers. The yield tensile stress generally decreases (from about 45 kPa to 5 kPa) as the fiber curvature increases (from 0 to 100) due to the reduction in fiber-fiber interlocking with smaller coordination numbers. Larger yield tensile stresses are observed for the samples with shorter length and smaller cross-sectional area, which involve fewer structural defects of weak fiber-fiber connections.
In the shear tests, the sample size effect also exists in that the smaller the sample length L s /d (<21) is, the larger the shear stresses are. The yield shear stress and coordination number both decrease with increasing fiber curvature κ for the S-shaped fibers, and the smaller yield shear stresses and coordination numbers are obtained for the U-shaped and Z-shaped fibers. The largest dimension of a fiber between two parallel tangential lines, namely the maximum Feret diameter, can be used to characterize the size of a fiber. It is observed that the yield shear stress exhibits an increase with the increase in the maximum Feret diameter for the fibers of various shapes. The yield shear stress can increase by an order as the normalized maximum Feret diameter increases from 0.55 to 1.
The results reported in this work may be useful for estimating how the mechanical responses of the fiber assemblies change when the original fiber shape (the steady shape under no external forces) changes, which can be applied to material designs with the fibers.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.