The mechanical behavior of knitted textiles is simulated using finite element analysis (FEA). Given the strong coupling between geometrical and physical aspects that affect the behavior of this type of engineering materials, there are several challenges associated with the development of computational tools capable of enabling physics-based predictions, while keeping the associated computational cost appropriate for use within design optimization processes. In this context, this paper investigates the relative contribution of a number of computational factors to both local and global mechanical behavior of knitted textiles. Specifically, different yarn-to-yarn interaction definitions in three-dimensional (3D) finite element models are compared to explore their relative influence on kinematic features of knitted textiles' mechanical behavior. The relative motion between yarns identified by direct numerical simulations (DNS) is then used to construct reduced order models (ROMs), which are shown to be computationally more efficient and providing comparable predictions of the mechanical performance of knitted textiles that include interfacial effects between yarns.

## Introduction

Knitted textiles are currently used as a new type of engineering material with manufactured internal architecture [13]. 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 [47], including shape predictions during the relaxation and draping processes [810], 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 [1114]. 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 [1518].

To improve computational efficiency, some nontraditional mesh-free techniques have also been reported [1922]. 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.

## Computational Approach

### Geometrical Model for Yarn Representation.

Many of the existing models that represent the shape of yarns in a knitted material [46] 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.

The challenge of producing physically accurate geometric models of yarns for knitted materials is framed in this approach as an optimization problem. Hence, a cost function $F()$ is defined, based on the shape of the yarn geometries, whose minimum represents a valid, physically realistic geometric initial condition for simulations. Function $F()$ consists of three terms defined in the below equation:
$F(Ci(Pj))=αFDistance+βFBending+γFLength$
(1)

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.

Fig. 1
Fig. 1
Close modal

The centerlines are specified with Catmull–Rom splines [29], which consist of C1 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.

The specific definitions of the three terms are as follows:
$FDistance=(D−2R)2$
(2)
where D is the shortest distance between the two splines that define the centerlines of two yarn models, R is the radius of the yarn, which also defines the offset surface of a spline tube, and the difference is squared in order to make the term continuous. In Eq. (2), D is approximated by computing numerous points on the yarn splines (in this case, 100 such points were chosen for each Bezier curve) and then finding the closest two points in these sets for the two separate splines. The remaining terms in Eq. (1) are defined as
$FBending=Fcont+β0∑kKjoin2$
(3)
where
$Fcont=∫‖C″(t)2dt$
(4)
and
$Kjoin={2σcos(θ2)π4<θ≤π−(π4)2ρθ+τ+π4ρ0≤θ≤π4$
(5)
In Eq. (4), $Fcont$ is the bending energy in the continuous interior portion of the individual cubic Bezier curves, defined as the integral of the square of the curve's curvature. In addition, $Kjoin$ approximates the curvature across the $k$ join points of the Catmull–Rom spline. See Breen et al. [10] for further details about this function and its parameters. In Eq. (3), $β0$ is a scalar that ensures that discrete bending energy $Kjoin2$ is proportional to the continuous bending energy $Fcont$. Now, $FLength$ is defined as
$FLength=∑l(Ll0−Ll)2$
(6)

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

### Interfacial Interactions in Direct Numerical Simulations.

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.

Fig. 2
Fig. 2
Close modal
In the relative movement between two yarn surfaces in each contact location under the tie constraint assumption, several nodes of the so-called slave surface are forced to have the same displacement as the corresponding node of the master surface based on the “surface-to-surface” discretization, shown in Fig. 2(a). In the case though of contact and friction, by taking the normal and tangent terms on the interfaces into consideration (Fig. 2(b)) and the principle of virtual work [30], one obtains
$δPC=δP+δG+δGT$
(7)
where $δPC$ and $δP$ are the virtual work of the whole system and of the deformable noncontact zones, while $δG$ and $δGT$ are the virtual work of the constraints in the normal and tangential direction of the contact surfaces, respectively. Hard contact in the normal direction is applied to prevent penetration and Coulomb friction is enforced in the tangent direction. These constraints are numerically implemented by a penalty method to improve computational efficiency, which rephrases the corresponding virtual work terms in the following way:
$δG=δGP=∫ΓcδγNpdΓ$
(8)
$δGT=∫TCδγT⋅tTdT=∫TCδγT⋅pμdT$
(9)

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.

When a relatively greater contact stiffness in Eq. (10) (where $N$ is the element shape function) compared to the element stiffness $KE$ in Eq. (11) is assigned, the displacement of the slave node approaches that of the master node and therefore the zero interpenetration can be approximated. The constraint equation enforced on slave nodes by this penalty method is expressed in Eq. (11) when no external force exists on the contact surfaces, where u and u* are the displacement of slave node and master node, respectively. Hence
$P=∫ΓcNpdΓ$
(10)
$(KE+P)u=Pu*$
(11)

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.

Table 1

Properties of contact and friction interactions

Contact-overclosure relationHard contact
Ratio of contact stiffness/element stiffness10
Friction modelIsotropic columbic friction
Ratio of stick stiffness/element stiffness10
Friction coefficient0.76 [32]
Contact-overclosure relationHard contact
Ratio of contact stiffness/element stiffness10
Friction modelIsotropic columbic friction
Ratio of stick stiffness/element stiffness10
Friction coefficient0.76 [32]

### Boundary Conditions.

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

Fig. 3
Fig. 3
Close modal

### Geometry Input for the Two-Dimensional Reduced Order Model.

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 xy plane can be obtained as shown in Fig. 4.

Fig. 4
Fig. 4
Close modal

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.

### Geometry Input for the Three-Dimensional Reduced Order Model.

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.

Fig. 5
Fig. 5
Close modal

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.

### Interfacial Interactions in Three-Dimensional Reduced Order Model.

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.

Fig. 6
Fig. 6
Close modal
Table 2

Kinematic constraints and properties used in the 3D-ROM

Interaction behaviorConstraints and properties
Tie constraints$uT=0,uR=0$
Isotropic friction(Course tension)
$kx=84.5 KN/m,ky=67.6 KN/m$, $kz=952.1 KN/m$
(Wale tension)
$kx=73.4 KN/m,ky=233.2 KN/m$, $kz=600 KN/m$
Interaction behaviorConstraints and properties
Tie constraints$uT=0,uR=0$
Isotropic friction(Course tension)
$kx=84.5 KN/m,ky=67.6 KN/m$, $kz=952.1 KN/m$
(Wale tension)
$kx=73.4 KN/m,ky=233.2 KN/m$, $kz=600 KN/m$

## Results and Discussion

### Local Interfacial Effects on Global Behavior.

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.

Fig. 7
Fig. 7
Close modal

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.

Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

### Geometrical Interpretation of Interfacial Effects on Global Mechanical Behavior.

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.

Fig. 10
Fig. 10
Close modal

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

Fig. 11
Fig. 11
Close modal

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.

### Computational Cost and Comparison With Beam Contact for Direct Numerical Simulations.

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 × 106, while the corresponding number of DOF of a 30 × 30 model dramatically rises to ∼20 × 106 for tie constraint and ∼53 × 106 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.

Fig. 12
Fig. 12
Close modal
Table 3

Degree-of-freedom and memory needed for different domain sizes (contact interaction)

c3w3c6w6c8w8c20w20c30w30
DOF311,9131.73 × 1063.30 × 1062.31 × 1075.31 × 107
Memory1158478810,39867,646154,825
c3w3c6w6c8w8c20w20c30w30
DOF311,9131.73 × 1063.30 × 1062.31 × 1075.31 × 107
Memory1158478810,39867,646154,825

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.

### Two-Dimensional Reduced Order Model Comparison With Direct Numerical Simulations.

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.

Fig. 13
Fig. 13
Close modal

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.

Fig. 14
Fig. 14
Close modal
Table 4

Line orientations of the geometry-based and MPS-based honeycomb with respect to the global x-, y-, z-axes

$1$$2$$3$$4$$5$$6$
2D ROM (honeycomb)(0, 90, 90)(66, 24, 90)(66, 24, 90)(80, 10, 90)(80, 10, 90)(0, 90, 90)
MPS-based (course)(1, 90, 89)(31, 59, 90)(36, 54, 90)(32, 58, 89)(33, 57, 89)(0, 90, 90)
MPS-based (wale)(1, 90, 89)(82, 8, 90)(82, 8, 90)(82, 8, 89)(82, 8, 89)(0, 90, 90)
$1$$2$$3$$4$$5$$6$
2D ROM (honeycomb)(0, 90, 90)(66, 24, 90)(66, 24, 90)(80, 10, 90)(80, 10, 90)(0, 90, 90)
MPS-based (course)(1, 90, 89)(31, 59, 90)(36, 54, 90)(32, 58, 89)(33, 57, 89)(0, 90, 90)
MPS-based (wale)(1, 90, 89)(82, 8, 90)(82, 8, 90)(82, 8, 89)(82, 8, 89)(0, 90, 90)

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.

### Verification of Reduced Order Model Use in Modeling the Mechanical Behavior.

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.

Fig. 15
Fig. 15
Close modal

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.

Table 5

Comparison of computational cost for a 3 × 3 basic knit with tie constraint

Element typeMesh size (m)Degrees-of-freedomMinimum memory required
3D full modelC3D8R3 × 10−4189,513128
2D ROMB315 × 10−443217
3D ROMB315 × 10−4599419
Element typeMesh size (m)Degrees-of-freedomMinimum memory required
3D full modelC3D8R3 × 10−4189,513128
2D ROMB315 × 10−443217
3D ROMB315 × 10−4599419

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.

Fig. 16
Fig. 16
Close modal

## Conclusion

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.

## Acknowledgment

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.

## Funding Data

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

## References

1.
Kiekens
,
P.
,
2012
,
Intelligent Textiles and Clothing for Ballistic and NBC Protection: Technology at the Cutting Edge
,
Springer
,
New York
.
2.
Struszczyk
,
M. H.
,
2014
, “
Design Aspects of Fibrous Implantable Medical Devices
,” International Congress on Healthcare and Textiles
, Izmir, Turkey, Sept. 25–26, pp. 84–92.
3.
Lorussi
,
F.
,
Galatolo
,
S.
,
Bartalesi
,
R.
, and
De Rossi
,
D.
,
2013
, “
Modeling and Characterization of Extensible Wearable Textile-Based Electrogonimeters
,”
IEEE Sens. J.
,
13
(
1
), pp. 217–228.
4.
Peirce
,
F. T.
,
1947
, “
Geometrical Principles Applicable to the Design of Functional Fabrics
,”
Textile Res. J.
,
17
(
3
), pp.
123
147
.
5.
Munden
,
D. L.
,
1959
, “
The Geometry and Dimensional Properties of Plain-Knit Fabrics
,”
J. Text. Inst. Trans.
,
50
(
7
), pp.
448
471
.
6.
Demiroz
,
A.
, and
Dias
,
T.
,
2000
, “
A Study of the Graphical Representation of Plain-Knitted Structures—Part II: Experimental Studies and Computer Generation of Plain-Knitted Structures
,”
J. Text. Inst.
,
91
(
4
), pp.
481
492
.
7.
Renkens
,
W.
, and
Kyosev
,
Y.
,
2010
, “
Geometry Modelling of Warp Knitted Fabrics With 3D Form
,”
Text. Res. J.
,
81
(
4
), pp.
437
443
.
8.
Kaldor
,
J. M.
,
James
,
S. L.
, and
Marschner
,
S.
,
2008
, “Simulating Knitted Cloth at the Yarn Level,”
ACM SIGGRAPH
2008, Los Angeles, CA, Aug. 11–15, p. 65.
9.
Yuksel
,
C.
,
Kaldor
,
J. M.
,
James
,
D. L.
, and
Marschner
,
S.
,
2012
, “
Stitch Meshes for Modeling Knitted Clothing With Yarn-Level Detail
,”
ACM Trans. Graph.
,
31
(
4
), pp.
1
12
.
10.
Breen
,
D. E.
,
House
,
D. H.
, and
Wozny
,
M. J.
,
1994
, “
A Particle-Based Model for Simulating the Draping Behavior of Woven Cloth
,”
Text. Res. J.
,
64
(
11
), pp.
663
685
.
11.
,
S. G.
,
Kallivretaki
,
A. E.
, and
Provatidis
,
C. G.
,
2007
, “
Mechanical Simulation of the Plain Weft Knitted Fabrics
,”
Int. J. Clothing Sci. Technol.
,
19
(
2
), pp.
109
130
.
12.
Cimilli
,
S.
,
Nergis
,
F. B. U.
, and
Candan
,
C.
,
2008
, “
Modeling of Heat Transfer Measurement Unit for Cotton Plain Knitted Fabric Using a Finite Element Method
,”
Text. Res. J.
,
78
(
1
), pp.
53
59
.
13.
Karimian
,
M.
,
Hasania
,
H.
, and
Ajeli
,
S.
,
2014
, “
Numerical Modeling of Bagging Behavior of Plain Weft Knitted Fabric Using Finite Element Method
,”
Indian J. Fibre Text. Res.
,
39
, pp.
204
208
.
14.
Rao
,
M.
,
Duan
,
Y.
,
Keefe
,
M.
,
Powers
,
B.
, and
Bogetti
,
T.
,
2009
, “
Modeling the Effects of Yarn Material Properties and Friction on the Ballistic Impact of a Plain-Weave Fabric
,”
Compos. Struct.
,
89
(
4
), pp.
556
566
.
15.
Carvelli
,
V.
,
Corazza
,
C.
, and
Poggi
,
C.
,
2008
, “
Mechanical Modelling of Monofilament Technical Textiles
,”
Comput. Mater. Sci.
,
42
(
4
), pp.
679
691
.
16.
Fillep
,
S.
,
Mergheim
,
J.
, and
Steinmann
,
P.
,
2017
, “
Towards an Efficient Two-Scale Approach to Model Technical Textiles
,”
Comput. Mech.
,
59
(
3
), pp.
385
401
.
17.
Yeoman
,
M. S.
,
Reddy
,
D.
,
Bowles
,
H. C.
,
Bezuidenhout
,
D.
,
Zilla
,
P.
, and
Franz
,
T.
,
2010
, “
A Constitutive Model for the Warp-Weft Coupled Non-Linear Behavior of Knitted Biomedical Textiles
,”
Biomaterials
,
31
(
32
), pp.
8484
8493
.
18.
Tanov
,
R.
, and
Tabiei
,
A.
,
2000
, “
Computationally Efficient Micromechanical Models for Woven Fabric Composite Elastic Moduli
,”
ASME J. Appl. Mech.
,
68
(
4
), pp.
553
560
.
19.
Wang
,
Y.
, and
Sun
,
X.
,
2001
, “
Digital-Element Simulation of Textile Processes
,”
Compos. Sci. Technol.
,
61
(
2
), pp.
311
319
.
20.
,
R.
, and
Robitaille
,
F.
,
2014
, “
Particle-Based Modeling of the Compaction of Fiber Yarns and Woven Textiles
,”
Text. Res. J.
,
84
(
11
), pp.
1159
1173
.
21.
Nocent
,
O.
,
Nourrit
,
J.-M.
, and
Rémion
,
Y.
,
2001
, “
Towards Mechanical Level of Detail for Knitwear Simulation
,”
J. WSCG
,
9
(1–3), p. 252.
22.
Cirio
,
G.
,
Lopez-Moreno
,
J.
, and
,
M. A.
,
2017
, “
Yarn-Level Cloth Simulation With Sliding Persistent Contacts
,”
IEEE Trans. Visualization Comput. Graph.
,
23
(
2
), pp.
1152
1162
.
23.
Vu
,
T. D.
,
Durville
,
D.
, and
Davies
,
P.
,
2015
, “
Finite Element Simulation of the Mechanical Behavior of Synthetic Braided Ropes and Validation on a Tensile Test
,”
Int. J. Solids Struct.
,
58
, pp.
106
116
.
24.
Leech
,
C. M.
,
2002
, “
The Modelling of Friction in Polymer Fibre Ropes
,”
Int. J. Mech. Sci.
,
44
(
3
), pp.
621
643
.
25.
Durville
,
D.
,
2012
, “
Contact-Friction Modeling Within Elastic Beam Assemblies: An Application to Knot Tightening
,”
Comput. Mech.
,
49
(
6
), pp.
687
707
.
26.
Chaouachi
,
F.
,
Rahali
,
Y.
, and
Ganghoffer
,
J.
,
2014
, “
A Micromechanical Model of Woven Structures Accounting for Yarn–Yarn Contact Based on Hertz Theory and Energy Minimization
,”
Composites, Part B
,
66
, pp.
368
380
.
27.
Bechtel
,
S. E.
,
Vohra
,
S.
, and
Jacob
,
K. I.
,
2012
, “
Contrasting the Predictions for Coulomb and Creep-Rate-Dependent Friction in the Modeling of Fiber-Draw Processes
,”
ASME J. Appl. Mech.
,
79
(
6
), p.
061001
.
28.
Liu
,
D.
,
Christe
,
D.
,
Shakibajahromi
,
B.
,
Knittel
,
C.
,
Castaneda
,
N.
,
Breen
,
D.
,
Dion
,
G.
, and
Kontsos
,
A.
,
2017
, “
On the Role of Material Architecture in the Mechanical Behavior of Knitted Textiles
,”
Int. J. Solids Struct.
,
109
, pp.
101
111
.
29.
Catmull
,
E.
, and
Rom
,
R.
,
1974
, “
A Class of Local Interpolating Splines
,”
Computer Aided Geometric Design
,
,
New York
.
30.
Belytschko
,
T.
,
Liu
,
W. K.
,
Moran
,
B.
, and
Elkhodary
,
K.
,
2013
,
Nonlinear Finite Elements for Continua and Structures
,
Wiley
, Hoboken, NJ.
31.
ABAQUS
,
2013
, “ABAQUS, Version 6.13, 2013 User's Manual,”
Dassault Systems
,
Pawtucket, RI
.
33.
Khandelwal
,
S.
,
Siegmund
,
T.
,
Cipra
,
R. J.
, and
Bolton
,
J. S.
,
2015
, “
Adaptive Mechanical Properties of Topologically Interlocking Material Systems
,”
Smart Mater. Struct.
,
24
(
4
), p.
045037
.