A Computational Approach to Model Interfacial Effects on the Mechanical Behavior of Knitted Textiles

Knitted textiles are currently used as a new type of engineering material with manufactured internal architecture [1–3]. Despite their widespread use and the existence of scalable production methods, the available tools for digital design coupled with performance predictions are relatively limited. In fact, although elaborate geometry models have been developed [4–7], including shape predictions during the relaxation and draping processes [8–10], the lack of physically accurate computational modeling related to the mechanical behavior of knitted textiles continues to pose challenges in engineering design using these new materials.

The complex internal structure of textiles, composed of geometrical unit cell arrangements, has resulted in simulation work focusing on their mechanical and multiphysics performance at the scale of the unit cells, primarily also due to the occurring computational cost when attempting to scale up such models [11–14]. In some investigations, multiscale analysis and computational homogenization were attempted to enable predictions, mostly focusing though on woven textiles because of their use as reinforcement phase in fiber-reinforced composites [15–18].

To improve computational efficiency, some nontraditional mesh-free techniques have also been reported [19–22]. Specifically, the concept of “digital element” was proposed to discretize yarns into pin-connected, digital-rod-element chains in which, though, physical properties were not assigned [19]. Furthermore, particle-based models were used to predict the compression behavior of textiles [20], while bounding volume models resulting from reduction of geometrical parameters were proven to be capable of simulating the dynamical process of a hanging fabric [21]. In addition, a force model with a sliding contact between yarns was demonstrated to result in significant reductions in computational cost [22]. However, most of these investigations were based on phenomenological analyses, with only weak or missing connection to actual three-dimensional (3D) textile behavior.

Meanwhile, the issue of interfacial effects at different scales in textiles is attracting great interest with the advent of new fibers and the use of advanced manufacturing methods. For example, a finite element analysis (FEA) approach to investigate the mechanical behavior of braded ropes involving contact and friction interactions between the internal yarns was recently developed [23]. Furthermore, particle-based models focusing on the compaction behavior between fibers have also been reported [20]. In addition, contact and friction between internal components of polymer fiber ropes were investigated to describe the rope deformations and loads [24]. Moreover, self-contact between beams was applied in a knot tightening simulation subjected to large deformations [25], while the equilibrium shape of a weave fabric accounting for yarn-to-yarn contact based on Hertz's theory under uniaxial and biaxial tension was investigated [26]. A creep rate-dependent friction model was also developed between filament tow and draw rollers to understand the behavior of fiber draw [27]. However, assumptions regarding the targeting geometry were made in these studies, which significantly decrease the complexity of the contact problem and, thus, their overall capabilities to realistically simulate the mechanical behavior of textiles.

In this context, this paper presents a novel computational approach for yarn-level interfacial effects in knitted textiles, which is shown to be capable of determining their influence on the global material response based on direct numerical simulations (DNS) using three-dimensional FEA models. Given the quantified interfacial forces extracted from these simulations, relations between the local geometry and the global response are also established. Based on these DNS results, a concept for reduced order modeling (ROM) is introduced to reduce the associated computational cost in DNS. Comparisons between the full FEA and ROM models are then performed, and the obtained results are evaluated based on their capability to predict the influence of local interfacial effects on the global mechanical behavior of knitted textiles.

Many of the existing models that represent the shape of yarns in a knitted material [4–6] are purely geometric in nature. Specifically, these models are based on defining parameterized control points that specify tube centerlines in 3D space. While these models may be sufficient for visualization purposes, in the sense that they present some possible geometric configurations of the constituent yarns in a knitted material, they are inadequate for physical simulations. Such models are further deficient because, being only geometric representations, they are devoid of any of the yarns' mechanical properties, and therefore in general define unphysical initial conditions for the subsequent simulations based on these geometries. These yarn models may, therefore, interpenetrate or bend in severe or unusual ways, making them unsuitable, e.g., for FEA simulations.

In order to produce valid initial geometric conditions for the simulations presented in this paper, an enhanced model of knitted yarns that incorporates mechanical properties with the underlying geometric representation of the yarns is implemented based on previous work by Liu et al. [28]. This approach produces initial geometric models that not only do not interpenetrate but also are in contact at traceable points while having a feasible and physically accurate overall shape.

where $F()$ is a function of the curves $Ci$ that define the centerlines of the tubes used to represent the yarns shown in Fig. 1.

The centerlines are specified with Catmull–Rom splines [29], which consist of C¹ piecewise cubic Bezier curves. The curves are defined by a set of control points $Pj$⁠, which are interpolated by the curves $Ci$⁠. Furthermore, $FDistance$ is defined in such a way that it is at a minimum when the distance between the spline tubes is zero, i.e., the yarn geometries are not intersecting, and are in contact at a point. $FBending$ defines the bending energy (as a function of curvature) of the yarn and ensures that the yarn geometry maintains a natural-looking shape. Finally, $FLength$ is the term that ensures that the yarn spline tube does not significantly change its length during optimization. The parameters $α, β, and γ$ are weights that allow for the modification of the strength of each term and provide the capability to adjust the importance and contribution of each component to the final result.

where $Ll0$ is the initial length of the lth Bezier curve in a yarn spline, computed as the distance between the first and last control points of the curve. In addition, $Ll$ is the current length of the curve, which is computed during the optimization as the control points move during this process. A subset of the control points $Pj$ is modified via an unconstrained quasi-Newton method in order to minimize $F()$⁠. The minimization in fact is constrained, not only by the choice of which control points are modified, but also by the fact that not all of the components of the modified control points are changed. By applying this approach in the textile unit cell in Fig. 1(a), which consists of three different yarns shown in red, green, and blue color, respectively, it is demonstrated that severe interpenetrations between yarns (represented by the white spheres at the contact locations between the blue yarn and green/red yarn in Fig. 1(a)) vanish in the optimized geometry shown in Fig. 1(b), which makes such digital representations more suitable for use within the FEA approach defined in Secs. 2.2–2.3.

Interactions between yarns in the DNS in this paper are assumed to include contact with friction, as opposed to tie constraints used previously by Liu et al. [28]. Figure 2 provides a schematic illustration of both the tie constraint (Fig. 2(a)) and the contact with friction (Fig. 2(b)) for comparison.

where $δγN$ is the virtual penetration rate in the normal direction and $δγT$ is the virtual relative velocity in the tangent direction. Furthermore, $p$ in Eqs. (8) and (9) stands for the contact pressure, $tT$ is the tangent traction, and $μ$ is the friction coefficient.

Coulomb friction is characterized by a stick-slide behavior. When the contact bodies are sticking, the tangential force is smaller than the critical frictional force and therefore similar constraints to Eq. (11) are applied to the slave nodes, with the contact stiffness $P$ being replaced by stick stiffness in this case. The critical and sliding frictions are related to the normal contact force by the friction coefficient, while the tangent tractions are always balanced along interfaces. For deformable objects in 3D space, tangential tractions can be represented by two orthogonal directions in the traction plane defined as $t1$ and $t2$ in Fig. 2(b). Numerically, an isotropic friction model assumes identical friction coefficients for the two tangential directions, while a so-called rough friction model does not allow any tangential movements. Since the yarns are the same everywhere and based on the fact that surface-to-surface discretization is adopted to smooth local concentrations [31], master and slave surfaces are chosen randomly in the full models.

For contact tracking, finite sliding is chosen to capture potentially severe interfacial interactions. In this case, the geometrical nonlinear analysis is activated and static stabilization is utilized to assist in convergence. The dissipated energy resulting from the stabilization is restrained to be less than 5% of the strain energy so that the accuracy of solution is not significantly impacted. Specific properties of the contact and friction parameters are listed in Table 1.

For the simulations in this paper, the yarn material was assumed to be isotropic polypropylene and a commercially available FEA code (abaqus 6.13/Standard) was used [28]. Based on the optimized geometry input in Fig. 1(b), in-plane stretches along the course (x-axis) and wale (y-axis) directions were performed by keeping one side of the domain fixed and applying prescribed displacements on the opposite side, while no constraints were applied on the other two sides of the domain. Linear hexahedron-reduced integration elements were used to avoid locking problems and achieve a higher computational efficiency. Additional details related to this FEA modeling setup can be found in Ref. [28].

Due to the geometry of the simulated knitted fabric and the free boundary conditions applied along the sides, the discrete yarns defined with contact and friction interaction on the top and bottom boundaries tend to slide when load is applied on the left and right boundaries under course tensile load. This effect prevents loading to be correctly transferred within the domain, while it additionally causes numerical convergence difficulties. Therefore, the top and bottom rows of yarns in the course tension simulations reported herein are defined by tie constraint in order to avoid any unnecessary rigid body motions, as shown in Fig. 3(a). For the same reason, the boundary conditions for the wale simulations are adjusted accordingly as shown in Fig. 3(b).

Based on the location of contact zones, e.g., in the FEA model of Fig. 3, a honeycomb structure can be obtained by directly linking the contact zones in the undeformed DNS model and by assuming that there are only small coordinate differences along the z-axis so that a 2D structure confined on the x–y plane can be obtained as shown in Fig. 4.

The 2D-ROM can then be discretized by linear beam elements with the same cross section as the yarns in the 3D DNS model. The joints formed in this approach at the intersecting points of the line segments in Fig. 4 represent the contact zones in which single nodes are shared by different elements. Essentially, this type of joint is equivalent to the tie constraint previously applied in the 3D DNS model, which restricts local kinematics. The effect that such simplification has on the mechanical behavior of knitted textiles is explored next.

To better capture the features of the full 3D DNS model, a 3D ROM is also constructed. Generally, a three-dimensional yarn can be represented by a beam in FEA, whose behavior depends on its centerline and cross-sectional area definitions. In the case of knitted textiles, though, the interfacial forces at contact locations have an important role in the overall mechanical behavior and therefore an approach is presented next to include such interactions in a 3D ROM. Specifically, to construct a 3D-ROM that contains information related to contact interactions, the following procedure was used based on the schematic in Fig. 5. First, contact points at the actual interface between two yarns were identified. Such contact points correspond to more than one node in the DNS model. For this reason, the contact points for the 3D-ROM described herein consist of points such as P1 and P2 in Fig. 5, which have coordinates computed by finding the centroid of the involved contact areas. In addition, a contact normal vector was defined in the direction of the resultant contact force on the yarn interface.

The contact points P1 and P2 in Fig. 5, in conjunction with the normal direction in each one of them, define the contact normal plane. It should be noted here that since this ROM definition is made in an undeformed configuration, the normal direction at points P1 and P2 is the same since loops are assumed symmetric. Furthermore, the intersection of the contact normal plane with the yarn centerlines was used to define points C1–C4 in Fig. 5, which consist of the pairs C1–C3 and C2–C4, each one of which contains a point along the centerline of the two yarns that relates to each contact points P1 and P2 of the 3D textile architecture. Such pairs along the centerlines of the yarns will be used to describe interfacial interactions in the 3D ROM models, as described next.

Based on the approach defined in Fig. 5, the 3D-ROM of Fig. 6 is obtained. To approximate the tie constraint, the relative motion between points C1–C3 and C2–C4 is eliminated. To simulate contact with isotropic friction, the interfacial interactions in the 3D-ROM are represented by equivalent linear springs. Hence, for a particular contact point in the DNS model (e.g., P1 or P2), the relative displacement between yarn interfaces and the net interfacial forces (contact forces along both the normal and tangential direction) are extracted, the ratio of which is used to define equivalent spring stiffness value for the x, y, and z directions, which are applied as forces between points C1–C3 and C2–C4. In addition, rotation at the contact points in the ROM is allowed since the rotational degrees-of-freedom (DOF) are inactive for solid element and consequently the rotational resistance is difficult to be obtained from the DNS model. The constraints and equivalent stiffness values used in the 3D-ROM are listed in Table 2.

The assumed yarn interaction effects defined in this paper are capable of capturing characteristics of the interfacial motion between yarns. To demonstrate this, Fig. 7(a) shows the computed evolution of the number of contact zones and their corresponding computed area estimation for one unit cell at the center of a model called a “3 × 3” in this paper since it consists of three unit cells in the course and wale directions as shown in Fig. 4. The computed global reaction force versus applied course strain is shown in Fig. 7(b) as a measure of global behavior, while the actual 3D FEA results that show the evolution of both contact zones and size are demonstrated in Fig. 7(c) in which the arrows are pointing to locations where actual contact is simulated to occur. The numbers 1–4 in Figs. 7(a) and 7(b) are associated with the four full field snapshots in Fig. 7(c). Note in Figs. 7(a) and 7(b) that slope changes in Fig. 7(b) occur right after point 3, which are associated with the increase of the number and size of the contact zones as shown in Figs. 7(a) and 7(c). Physically, this increase of the contact zones corresponds to the entanglement of the yarns as they are pulled in the course direction, as explicitly shown in Fig. 7(c). Specifically, no contact zones exist in the first snapshot of Fig. 7(c), which corresponds to a global course strain of 0.05%. This was expected as all interpenetrations were diminished by the optimization approach, previously described. As further course tensile load is applied, contact starts to occur in symmetric locations of the unit cell (see the third snapshot in Fig. 7(c)) and as a result, the number of contact zones is increasing to 4 (at 1.95% of applied strain). This contact topology remains constant until the total 5% strain is applied in this simulation.

A complex mechanical behavior under in-plane loading conditions was previously reported by Liu et al. [28], which motivated the approach presented in this paper. Specifically, it is shown in Fig. 8 that a pronounced out-of-plane motion, demonstrated by comparing the deformed shape with the original undeformed configuration in addition to plotting the transverse (y-axis/wale) contraction, occurs during course (x-axis) tension loading for the contact and friction model. Interestingly, the out-of-plane displacement is symmetric with respect to the free edges (top and bottom), while the interfacial forces at each contact point can be directly computed, as shown in Fig. 8. Hence, in accordance with the previous work by Liu et al. [28], the yarn interfaces are found to have an important role in the load transfer within the knit material architecture.

In Fig. 9, reaction forces (Fig. 9(a)) and out-of-plane displacement (Fig. 9(b)) are compared for the different assumptions of interfacial interactions under the same course tensile loading. Based on the results in Fig. 9(a), it could be argued that if the slope of the curve defined as reaction force versus applied course strain is a measure of overall stiffness of the knitted textile, then the tie constraint represents the stiffer case while the isotropic friction is the softer one. This result can be explained by the fact that in the cases where the normal and tangential components of friction vary, there are more DOFs for the yarn to move which results to additional ways for the applied load to be carried and transferred. Since energy is expended in this type of interaction, the material becomes structurally less stiff. Based on the results in Fig. 9 and similar results obtained but not shown in this paper, a smaller friction coefficient value does not impact the conclusion drawn, which relates to the fact that the stiffer behavior is demonstrated in the case of hard/rough friction.

The interfacial forces acting at all yarn contact locations, as shown in Fig. 8, are composed of the normal contact force and the tangent friction (shear) force based on the assumptions previously made regarding contact and friction. To further explore the effect of interfaces on local yarn kinematics, individual yarns in knitted textiles are considered as parts of an interconnected system of yarns. In this context, a free body diagram of each yarn based on the model shown in Fig. 8 could be formed as shown in Fig. 10(a). Based on this description, each yarn is subject to a net interfacial force exerted at all FEA nodes of each contact zone that can be resolved in components aligned with the Cartesian axes. Note that for the linear hexahedron element used in the study, no rotations were prescribed as nodal degrees-of-freedom. In addition, a reaction force as well as a moment is created at each point where the yarn is interlocked with the rest of the knitted textile structure.

Given the FEA results in Fig. 10(a) and by adopting a contact point of view to describe the local motion, the path schematically shown in Figs. 10(b) and 10(c) was observed, which includes a rotation between states 1 and 2 followed by a slide from state 2 to a new state 3. In addition to this path observed in each yarn, it was found that the contact areas also spin around their main axis as noted in Fig. 10(c). This is the first time in the authors' knowledge that these intricate local motion characteristics of yarns in knitted textiles are simulated and physically explained as shown in this section.

Based on the results in Fig. 10, simulation results are presented with varying course spacing $(C)$ in Fig. 11(a) and yarn diameter in Fig. 11(b), since these two geometrical properties are expected to affect the complex motion based on the motion description in Fig. 10. Indeed, the out-of-plane displacement is found to increase with increasing course spacing in Fig. 11(a), as in this case, the moment axis that causes the motion from state 1 to state 2 in Fig. 10(b) also increases. In addition, an increase in the yarn diameter increases its moment of inertia, which reduces the capability to spin and therefore the out-of-plane displacement decreases as shown in Fig. 11(b).

Consequently, the free body diagram description in Fig. 10 assists in interpreting quantitatively the local and global kinematics of knitted textiles and verifies that the modeling approach followed is sensitive enough to address practical design inputs to manufacturing as both the diameter of the yarn and the course spacing can be controlled by existing manufacturing equipment.

The DNS presented in this paper were based on models composed of 3D solid yarns, which assisted in the investigation of yarn interfacial properties. However, such full field numerical investigations are also computationally intensive, especially when the computational domain size is increasing. In Fig. 12(a), the total DOF of a 3 × 3 knit is found to be ∼0.2–0.3 × 10⁶, while the corresponding number of DOF of a 30 × 30 model dramatically rises to ∼20 × 10⁶ for tie constraint and ∼53 × 10⁶ for contact/friction assumptions. The corresponding memory needs increase as listed in Table 3. It should be noted that it was previously demonstrated by the authors that knitted textile simulations show significant size-dependent behavior when treated as free standing structures, which combined with the demonstrated here computational cost motivate herein the development of reduced order models.

Additionally, a 3 × 3 textile model employing beam elements and beam-to-beam contact was compared to the full DNS result before proceeding to reduced order model techniques. It is shown in Fig. 12(b) that the beam contact model demonstrates a softer behavior compared to full model as the cross-sectional deformation is ignored and only one contact node is captured for nonparallel beams. The fact that beam-to-beam contact has an effect on textile behavior further contributes to the demand for a reduced order modeling technique.

To check the validity of 2D (honeycomb) ROM constructed by simply connecting contact zones as shown in Fig. 4, volume averages of stresses for each section in Fig. 13 were used and the maximum principal stress (MPS) values for each section was calculated similar to the so-called thrust-lines method used previously for topologically interlocked materials to obtain ROM definitions [33]. For both course and wale tension cases, Fig. 13 shows that the computed direction of the MPS at different locations varies; however, it can be observed that it tends to align with the direction of the yarn loops. Intricate flows with complicated directions appear around the contact areas, confirming that contact points direct, transfer, and can deflect the stress flow in knitted textiles.

Given the computed directions of the MPS shown in Table 4, two additional honeycomb structures (named MPS-based course/wale) were constructed by setting the same distance between contact zones along the x-axis and y-axis and using a set of three angles to position each of the six segments in space as shown in Fig. 14. Based on the data in Table 4 and the schematics in Fig. 14, differences were observed between the ROM obtained based on geometry in Fig. 4 and MPS calculations. Specifically, it was found that the ROM based on the calculated MPS lines is compressed for course tension loading and stretched for wale tension loading along the wale direction. This behavior is caused by the tendency of the MPS lines to align with the externally applied load creating an issue of which of the two possible ROM representations should be chosen to represent any arbitrary loading. Therefore, the ROM structure obtained using the geometry-based approach seems to be a trade-off between physical accuracy and numerical implementation.

This choice is actually and in part also justified by the data in Table 4 from which it can be concluded that the ligaments of the produced honeycomb structure defined based on geometrical considerations are consistent with corresponding directions formally computed by the concept of MPS lines.

Comparisons of the reaction force versus applied strain curves for both the course and wale directions of the ROM and 3D full FEA models are shown in Fig. 15. It is clearly seen in Fig. 15(a) that the 3D-ROM approaches much closer the behavior of the DNS for both tie and isotropic friction definitions under both course and wale direction tensile loading.

Regarding computational costs, for a 3 × 3 model, the total number of DOF decreases from 189,513 to 5994 and the memory required is reduced to 19 MB compared to 128 MB for the DNS results, as shown in Table 5. This is a crucial parameter to be considered as Fig. 12 shows the computed increase of the DOF and required memory for DNS as a function of computational size in the case of contact and friction, which justifies the need for ROM in the case of knitted textiles. In addition, it should be noted that although obviously there is an increased number of DOF in the case of 3D-ROM compared to the 2D-ROM, the corresponding increase in terms of memory requirements is minimal, predominantly because of the use of beam elements. Similar computational cost benefits were seen in the case of isotropic friction.

Regarding the out-of-plane motion observed for example in Fig. 8 using the DNS approach, the results shown in Fig. 16 demonstrate that the localizations computed in a 3 × 3 model using DNS are actually captured also by the 3D-ROM. In fact, this out-of-plane motion is lost when using the 2D honeycomb ROM due to its in-plane nature. Furthermore, it is also observed in the 3D-ROM results that the maximum out-of-plane motion in the case of isotropic friction is larger with respect to the tie constraint case, which agrees with the corresponding trend observed in DNS results.

Numerical simulations were performed on knitted textiles based on geometric representations in which yarn orientation were optimized by a minimization process. The computational results obtained can explain the evolutionary interfacial interactions during tensile loading. Specifically, the interactions between yarns were investigated by adopting a contact definition with friction, which introduces a more physically realistic interfacial activity compared to previously reported tie constraints. The contact with friction definition was found to produce contact forces along the contact normal direction, as well as frictional forces in the corresponding tangent plane, which resulted in a physically accurate decomposition of the type of motion performed by yarns in knitted textiles. Such a direct numerical approach was found, however, to be computationally expensive with increasing domain size and therefore a process to define both two- and three-dimensional reduced order computational models was presented. The comparison between direct and reduced order results showed good agreement, which can be exploited in the development of numerical tools for physics-based modeling of knitted textile behavior.

The authors thank Professor Vadim Shapiro from the University of Wisconsin-Madison for useful discussions related to the research presented in this paper. The numerical results reported were obtained by using computational resources supported by Drexel's University Research Computing Facility. In addition, assistance in obtaining some of the computational results and interpreting them was provided by Dr. Seid Koric from the National Center for Supercomputing Applications (NCSA) at the University of Illinois at Urbana-Champaign.

• Division of Civil, Mechanical and Manufacturing Innovation (Grant Nos. CMMI 1344205 and CMMI 1537720).