Programming Mechanics in Knitted Materials, Stitch by Stitch

We performed experiments on eight types of yarn that are classified in three categories: (1) two are large-gauge yarns (9-12 wraps per inch, WPI), (2) five are fine-gauge yarn (30-40 WPI) , and (3) the yarn used for the therapeutic glove prototype (14-18 WPI) (see the Knitted Glove Prototype section).

We used Brava worsted yarn (28455-White) from KnitPicks™, which is 100% acrylic yarn, hereafter referred to as the “acrylic yarn” and 082L Pearl cotton 3/2 (color 1800-13 sapphire) from Halcyon Yarn™, which is 100% cotton yarn, hereafter referred to as the “cotton yarn.” For each of the types of fabrics, we recorded the average yarn diameter within the fabric stitches as well as the average yarn lengths per stitch. Supplementary Tables 2 and 3 display the measurements for the acrylic and cotton yarn, respectively. We measured the bending rigidity, an approximate interaction potential, and the stress versus strain relationship for both types of yarn. We perform four uniaxial experiment runs on the samples to obtain the stress versus strain relationship.

We used a Taitexma™ Industrial Knitting Machine to create four types of fabrics with both the acrylic and cotton yarn: stockinette, garter, $1\times 1$ rib, and seed. Each fabric sample consisted of 31 rows and columns and were made at equal tensions and stitch size settings on the machine (Supplementary Note 11).

For an accurate model development, we obtained finer details of the fabric stitches. We created smaller copies of the experimental samples. We used a caliper to measure the average diameter of the yarn in situ. We then dissected them to obtain average yarn lengths per stitch for the four types of fabric.

We additionally fabricated five sets of samples (where each set contained the four types of fabrics) made from different lace weight yarns. Of the five, three yarns were from ColourMart™: heavy lace weight alpaca mohair silk mokka 811 ecru, heavy lace weight kid mohair and silk special celeste, and 2/28NM lace weight cashmere 8l brume (beige) each referred to as “lace-weight alpaca mohair”, “lace-weight blue mohair”, and “lace-weight cashmere” respectively. The other two lace weights were Bambu 12 Gauge 100% Bamboo in the color 010 Rice from Silk City Fibers™, hereafter referred to as “lace-weight bamboo”, and Tamm Petit 2/30 T4201 White 100% acrylic yarn from The Knit Knack Shop™, hereafter referred to as “lace-weight acrylic.” These samples were fabricated on a STOLL CMS 530 HP Industrial Knitting Machine and each contained 32 rows and 32 columns. Stockinette and garter were made with a stitch size setting of 12 while rib and seed were made at size 11 (Supplementary Note 11 and Supplementary Fig. 14). All other machine parameters were kept the same. Each sample was fabricated twice with buffer regions either along its vertical or horizontal axis to aid with the uniaxial stretching experiments.

Similar to the acrylic and cotton yarns mentioned above, we measured the bending rigidity for each of these yarns and extract the stress versus strain relationship via uniaxial experiments. Five experiment runs were performed on each sample.

To obtain the length of yarn per stitch for the lace weight samples, each sample was weighed and, using the mass density for the yarn types, the average length of yarn per stitch was estimated.

To perform the uniaxial stretching experiments, we designed a setup such that fabric samples had external forces uniformly applied to the boundary. We 3D printed clamps to use on both ends of the fabric samples and then had a dynamometer hooked on to one of the clamps that could be moved with a threaded rod. All components of the experiment were designed to move on guiding rails to keep everything level and prevent lateral and torsional motion. We designed the clamps with several teeth to effectively hold down both ends of the fabric sample and prevent slipping.

For each sample, we clamped the fabric on opposite ends. During the experiments, we positioned and leveled a camera above the sample. Colored pins were placed in the fabric and red points were painted on the clamps to aid with tracking during the analysis. The dynamometer was zeroed before the experiment and then incrementally moved by turning the threaded rod, applying the external force $F_{x}$ (or $F_{y}$) to the sample boundary until reaching its maximum force (30 N). Experiments were performed slowly, to approximate a quasistatic regime, stretching from a relaxed configuration to maximum extension over 1-3 minutes. An initialization is done for each sample where they are run through the entire experiment. This run is not included in the presented data as it is meant to break apart initial fiber connections and handling bias. Between subsequent experiment runs, the fabrics were reset to their initial resting length and briefly stretched in their transverse direction. We then waited five minutes before the next experiment run. We performed experiments along both axes of the fabrics (along its $x$- and $y$-direction).

We looked at the uniaxial response by tracking the length and waist dimensions as the external force is exerted on the boundary. For the overall bulk response, we used Fiji (https://imagej.net/Fiji) image processing software with the TrackMate plugin to track the pins and clamps on each of the sample videos and analyzed the position change of the coordinates (see Supplementary Fig. 1). The dynamometer reading was recorded using optical character recognition (OCR) on its seven segment display, ensuring stress and strain data were synchronized. Raw experimental data for the experimental runs on the acrylic, cotton, and therapeutic glove fabric samples can be found in Supplementary Fig. 3.

For the uniaxial experiments on the lace weight samples and the therapeutic glove swatches, we use an Instron Universal Testing Machine (UTM) Model 68SC-1. We 3D printed unique clamps with teeth to fit into the machine grips and ensure no slip boundary conditions during testing. A camera is focused on the sample while stretching and two pins are placed along the transverse direction. The clamp separation is measured and the displacement is tracked by the Instron software. The samples are stretched at a rate of 0.5 mm/sec and the lace weights are stretched to 25 N while the glove samples are stretched to 30 N. The remaining experimental procedure and analysis is the same as detailed previously except that the force data was synced by matching time steps from the tracked transverse data (acquired with Fiji) and the force data (acquired with the UTM).

For one uniaxial experiment, we captured the nonaffine displacement fields throughout the entire sample under stress. We clamped the acrylic garter sample and dusted graphite powder to create a speckle pattern for tracking aid. The camera was again leveled above the sample but positioned closer to capture more detailed deformation. To analyze and track the displacement fields we use the 2D digital image correlation (DIC) MATLAB software, Ncorr (https://www.ncorr.com/).

Yarn has a hierarchical filamentous structure with internal stresses and friction arising from the manufacturing process that complicates determining a bending modulus through cantilever experiments. Since probe-based measurements, such as the three-point flexural test, inevitably lead to compression of the yarn’s cross-section, we find that cantilever experiments yield the most consistent results, using simple approximations to the yarn shape. A schematic of the setup is shown in Supplementary Fig. 4.

Looking at four increments of yarn length ranging from 10 cm to 25 cm, we cut out five samples at each length and perform bending experiments on the yarn. For the lace weight yarns, the lengths cut were 6 cm, 9 cm, 10.5 cm, 12 cm, and 15 cm. Each sample is cut with an additional 10 cm of yarn that is adhered on a flat surface with double-sided tape. The yarn is hung off the edge of the platform and bends under its own weight due to gravity, adopting an approximately parabolic shape. A camera is positioned level to the setup and images the yarns’ behavior. We apply a blur and binarize filter to the images to isolate the yarn. Taking the points that compose the yarn shape, we fit a 4^th degree polynomial to find the approximate centerline of the yarn (see Supplementary Note 2).

We used a Zwick/Roell Z010 Universal Testing Machine (UTM) to perform compression experiments on the yarn. Three yarn samples of length 20 mm (for the acrylic yarn) and 30 mm (for the cotton yarn) were compressed between a probe tip of 5 mm in diameter and a custom acrylic stage also 5 mm in diameter. The UTM probe tip was slowly lowered onto the yarn, resulting in quasistatic measurements of the restoring force as a function of probe height. A schematic of the setup is shown in Supplementary Fig. 5.

To simulate the equilibrium configurations of knitted stitches, we modeled yarn as inextensible elastica with bending modulus $B$ and fixed total length $L$ per stitch. Interactions between overlapping yarn were treated with a hard-core, soft-shell model with a functional form derived from experimental measurements. Equilibrium configurations were determined by numerically minimizing the total yarn energy, given by the sum of the bending energy $E_{\rm bend}=(B/2)\int_{0}^{L}{\rm d}s\,\left|\partial_{s}\hat{\mathbf{t}}\right|^{2}$ and the core-shell interaction energy $V_{\rm int}$, with a fixed total length constraint. To perform this numerical minimization, we represented yarn configurations as degree-5 Bézier spline curves with degrees of freedom encoded by a collection of Bézier curve control points. The resulting space curves are twice continuously differentiable with respect to its arclength parameter $s$, and thus have continuous curvature. For more details, see Supplementary Note 4.

To craft the knitted glove prototype, we used Rowan™ Baby Cashsoft Merino which is composed of 57% wool, 33% acrylic, and 10% cashmere. We used four different colors to knit the four types of fabric in the glove: blue for stockinette, orange for garter, green for rib, and pink for seed.

The fabric was knit by hand (see Supplementary Note 11) on US size 2 needles (2.75 mm in diameter) except for the stockinette regions which were knit on US size 0 needles (2 mm in diameter).

Miniature test swatches were made of each type of fabric to assist with glove design and to perform uniaxial stretching experiments on. A stockinette sample was knitted on US size 2 needles for comparison. Each sample underwent five experimental runs on the UTM. Each sample has 25 columns and 34 rows.

From the uniaxial stretching experiments, we extracted the average fabric response by expressing the applied force $F_{x}$ (or $F_{y}$) in terms of stress components $\sigma_{xx}=F_{x}/W_{y}$ (or $\sigma_{yy}=F_{y}/W_{x}$), where $W_{x}$ and $W_{y}$ are fabric widths measured at the clamped edges.

To obtain the average strain response of the fabric, we focused on the displacement of four points placed on the axes of symmetry of the fabric. Two points were marked with pins located on the transverse axis and the other two points were marked by painted dots on the clamps along the direction of stretching. In the case where the fabric is stretched along the $y$-direction, two red points painted on the clamps are aligned such that the line connecting them lies along the y-axis and has length $L_{y}$. The other two pins were positioned close to the waist of the fabric such that the line connecting them lies along the $x$-axis and has length $L_{x}$. For the lace weight and the glove prototype samples, there are only two pins along the transverse axis. The clamp displacement is tracked via the UTM. The principal components of the strain tensor are then $\varepsilon_{xx}=(L_{x}-L_{x,0})/L_{x,0}$ and $\varepsilon_{yy}=(L_{y}-L_{y,0})/L_{y,0}$, where $L_{y,0}$ and $L_{y,0}$ are the respective pin separations of the un-stretched fabric (see Supplementary Fig. 5).

Figure Supplementary Fig. 6 shows the results of the stress-strain analysis for the cotton yarn.

The yarn of linear mass density $\lambda$ and length $L$ is adhered to the edge of a flat surface at $x=0$ and $y=0$ and the free end is allowed to drape under gravity. Equilibrium distributions of the internal moment $\mathbf{M}(s)$ and force $\mathbf{T}(s)$ are governed by the Kirchhoff rod equations LandauLifshitz1986 ,

	$\displaystyle\partial_{s}\mathbf{T}(s)-\lambda g\hat{\mathbf{y}}$	$\displaystyle=0\,,$		(Supplementary Equation 1a)
	$\displaystyle\partial_{s}\mathbf{M}(s)+\hat{\mathbf{t}}(s)\times\mathbf{T}(s)$	$\displaystyle=0\,,$		(Supplementary Equation 1b)

where $s\in[0,L]$ is the arclength coordinate with $s=0$ at the fixed end and $s=L$ at the free end. Integrating Eqs. Supplementary Equation 1a, Supplementary Equation 1b and using the free end boundary conditions $\mathbf{T}(L)=\mathbf{M}(L)=0$, we have

	$\displaystyle\mathbf{T}(s)$	$\displaystyle=\lambda g(s-L)\hat{\mathbf{y}}\,,$		(Supplementary Equation 2a)
	$\displaystyle\mathbf{M}(s)$	$\displaystyle=\lambda g\hat{\mathbf{y}}\cdot\int_{s}^{L}{\rm d}s^{\prime}\,\hat{\mathbf{t}}(s^{\prime})(s^{\prime}-L)\,.$		(Supplementary Equation 2b)

Evaluating Supplementary Equation 2b at the clamping point $s=0$ and integrating by parts, we find

$$\mathbf{M}(0)=-\lambda Lgx^{*}\hat{\mathbf{y}}\,,$$

(Supplementary Equation 3)

where $x^{*}=L^{-1}\int_{0}^{L}{\rm d}s\,x(s)$ is the $x$-component of the center of mass of the hanging yarn, recovering the basic result that the total gravitational torque applied to the yarn at the fixed boundary is simply the total gravitational force of the yarn, $-\lambda Lg\hat{\mathbf{y}}$, times its lever arm $x^{*}$. Finally, assuming the linear constitutive relationship $M(s)=B\kappa(s)$, we solve for the bending modulus $B$ of the yarn, viz.

$$B=\frac{\lambda Lgx^{*}}{\kappa(0)}\,,$$

(Supplementary Equation 4)

where $\kappa(0)$ is the curvature discontinuity of the yarn at the clamping point. We perform experiments on yarn samples of varying length (see Supplementary Fig. 8).

In order to extract the curvature discontinuity at the suspension point, we must determine the shape of the yarn in space. After an image is taken and cropped to the suspended yarn, the exposure is adjusted to maximize contrast between the yarn and the background. We apply a blur and binarize filter to the yarn images and fit the white pixels in the domain $y>-y_{\rm max}/5$ to the 4th-order polynomial curve $y(x)=(a/2)x^{2}+(b/4)x^{4}$. This is reminiscent of how Cornelissen and Akkerman corn used a polynomial fit to study yarn deflection during cantilever experiments. Images are blurred based on how many pixels out we can see stray fibers. For yarns containing certain fibers, stray filaments will inhomogeneously extend many yarn-radii from the spun center of the yarn. For these exceptionally fuzzy yarn types (alpaca mohair and blue mohair), the blur required to erase stray fibers is so great that a reliable fit to the core of the yarn is not feasible. In these cases, a simple blur filter is not sufficient in isolating the core and, after adjusting the exposure, the outer halo of the yarn was manually painted out before filtering the image. Using the polynomial fit, we extract the curvature discontinuity $\kappa(0)=a/(1+a^{2})^{3/2}$. The $x$-component of the center of mass, $x^{*}$, is approximated by the average over all $x$-coordinates of the binarized image of the yarn.

The results of these fits for all yarn types used in this study are reported in Supplementary Table 1.

Using force vs. probe height compression measurements from the yarn compression experiments (see Supplementary Fig. 9), we find an effective yarn interaction potential energy for use in the simulations. The yarn compression data show nonlinear behavior of the yarn’s resisting force as a function of probe height, exhibiting a soft regime for low compression that stiffens as the constituent fibers are forced to pack into a small volume for high compression. Noting that the stress versus strain measurements of the fabric attain maximum stresses of approximately 1 N/mm, we argue that it is sufficient to find an approximate force versus compression depth that follows the yarn compression data up to 1 N. We fit the compression data to a model force law given by

$$F_{\rm probe}(\rho)=A_{\rm comp}k\frac{\rho_{0}}{p}\left[\left(\frac{\rho_{0}}{\rho}\right)^{p}-1\right]$$

(Supplementary Equation 5)

for forces between $10^{-3}$ N and 1 N, where the lower bound was chosen to cut out fluctuations on the measured force presumably due to the corona of wispy fibers sticking out of the yarn. Note that this form is similar to the contact interaction assumed by Kaldor 2008 kaldor8 , except that the exponent $p$ is left as a fitting parameter. Here, $\rho=\delta_{\rm c}-\delta$ is the thickness of the yarn when the probe is at depth $\delta$, where $\delta_{\rm c}$ is a cutoff depth, modeling an effective incompressible “core” of the yarn. We use the probe depth at 3 N as the cutoff height. This model was chosen because it captures the compression-stiffening behavior of the yarn for low to moderate compression, where $p>1$ is a fitting parameter that encodes this nonlinear behavior. Furthermore, $F_{\rm probe}\rightarrow 0$ as $\rho\rightarrow\rho_{0}$, where $\rho_{0}=\delta_{c}-\delta_{0}$ is the uncompressed thickness of the yarn, where $\delta_{0}$ is the probe depth at the edge of the yarn (here taken to be when $F_{\rm probe}\approx 10^{-3}$ N). The fitting parameter $k$ sets the scale of the yarn’s compressional rigidity per area, with $A_{\rm comp}$ representing the compressed area of the yarn, which we approximate as the diameter of the UTM tool (5 mm) times the diameter of the yarn being compressed. For small compressions, where $\rho-\rho_{0}=\Delta\rho\ll\rho_{0}$, the compression force is approximately $F_{\rm probe}\simeq-A_{\rm comp}k\Delta\rho$ and $A_{\rm comp}k$ measures an effective spring constant. The results of these fits are shown in Supplementary Table 4.

In the simulations, two yarn segments in contact mutually compress each other. We approximate the compression in terms of the centerline-to-centerline distance of two yarn segments, $R(s,s^{\prime})\equiv|\bm{\gamma}(s)-\bm{\gamma}(s^{\prime})|$, where $\bm{\gamma}(s)$ and $\bm{\gamma}(s^{\prime})$ are two centerline points. For fixed centerline points, the compressed thickness is taken to be $\rho=R-2r_{\rm core}$, where $r_{\rm core}$ is an effective “core radius” of the yarn, representing the hard core cutoff radius, and $\rho_{0}=2r-2r_{\rm core}$, where $r$ is the outer radius of the yarn. This type of soft-shell, hard-core model has been used previously in the simulation method of Sperl 2022 sperl_estimation_2022 . The interaction potential energy density is given by

$$\mathcal{V}_{\rm int}(\zeta)=\left\{\begin{array}[]{cc}k\frac{(2r-2r_{\rm core})^{2}}{p(p-1)}\left[\zeta^{1-p}-1-(p-1)(1-\zeta)\right]&\,\,\text{for}\,\,\zeta<1\\ 0&\,\,\text{for}\,\,\zeta\geq 1\end{array}\right.\,,$$

(Supplementary Equation 6)

where $\zeta\equiv(R-2r_{\rm core})/(2r-2r_{\rm core})$ is a non-dimensional measure of compressed thickness.

For a more physically accurate model of compression in the future, we would like to study compression in the method that Park and Oh park_bending_2006 developed for bending, which takes into account the hierarchical structure of spun yarn.

In the uniaxial stretching experiments, the fabrics attain an hourglass-like waist as the force acting along one direction causes a response in the transverse direction. This response is a result of the fabric’s Poisson effect, i.e. stretching the fabric in one direction thins it in transverse directions. For example, if the fabric is stretched in the $x$-direction so that the strain component $\varepsilon_{xx}$ is positive, then it thins in the $y$-direction. Due to the constraints on transverse deformation imposed by clamps on two of the edges, the strain field component $\varepsilon_{yy}<0$ varies along the fabric, and the magnitude of the strain reaches a maximum where the waist narrows. The pin tracking approach we employed measures the deformation in this region.

Due to the anisotropy of each fabric, the transverse deformation is characterized by two Poisson ratios, $\nu_{yx}$ and $\nu_{xy}$. In principle, these ratios describe the linear response of the fabric under two different deformation protocols: if the fabric is stretched in the $x$-direction (i.e. $\varepsilon_{xx}>0$) then the transverse response is $\varepsilon_{yy}=-\nu_{yx}\varepsilon_{xx}$; if the fabric is stretched in the $y$-direction then the transverse response is $\varepsilon_{xx}=-\nu_{xy}\varepsilon_{yy}$. We measured these ratios for experimental and simulation results, and they are reported in Supplementary Table 10, Supplementary Table 12 and Supplementary Table 20.

Here, we provide a scaling rationale for the constitutive relationship. The results presented here are consistent with the RS model of SUPPLEMENTARY NOTE 8, including the form of the strain-stiffening term. First, consider the various lengthscales that describe a knit stitch. These include the yarn radius $r$, the length of yarn per stitch $L$, and a “bending lengthscale” $\ell_{\rm b}\sim\sqrt{B/T}$, where $B$ is the yarn’s bending modulus and $T$ is the tension of the yarn, which is obtained from dimensional analysis. In fact, $\ell_{\rm b}$ has a simple physical interpretation: if one considers an arc of radius $R$, the work done in stretching the arc’s radius to $R+\delta R$ under constant tension scales as $\sim T\delta R$, which is counteracted by a change in the bending energy, which scales as $\sim-B\delta R/R^{2}$, and the two generalized forces are in equilibrium if $R\sim\ell_{\rm b}\sim\sqrt{B/T}$. Thus, $\ell_{\rm b}$ can be regarded as the radius of curvature that dominates the bending energy of a curve in mechanical equilibrium under tension. As $T$ increases, $\ell_{\rm b}$ decreases, so that the bending energy of a stitch is increasingly concentrated to small regions of high curvature, which must be located at the entangled regions. Note that in the arc approximation, the $\ell_{\rm b}$ is both the radius of curvature and the arclength of the curved regions, so that $E_{\rm bend}\sim 1/\ell_{\rm b}$. Since the radius $r$ of the yarn represents a lower limit of the radius of curvature and the stitch length $L$ determines the periodicity of yarn shape, $\ell_{\rm b}$ is bounded by these two lengths: $r\lesssim\ell_{\rm b}\lesssim L$.

The limits of this region correspond to two physical cases. Case one, $\ell_{\rm b}\approx L$, occurs when there is little to no applied external force. Under small strains $\varepsilon$, the bending length decreases linearly, such that changes in $\ell_{\rm b}$ go as $\Delta\ell_{\rm b}\sim-\varepsilon\ell_{{\rm b}_{0}}$, where $\ell_{{\rm b}_{0}}$ is the bending length at zero strain. This results in the linear stress-strain relationship $\sigma_{\rm low}\propto\varepsilon$. Case two, $\ell_{\rm b}\approx r$, occurs when the fabric is under high external load. Here, the yarn segments within each crossover region clasp increasingly tightly around each other. In this regime, the total length of yarn $L$ can be approximated as $L\approx c(\lambda_{\rm max}+r)$, where $c$ is a numerical prefactor of $\mathcal{O}(1)$ and $\lambda_{\rm max}\approx(1+\varepsilon)\lambda_{0}$, the maximal separation between crossover regions, varies linearly with the average separation between crossover regions in the unstrained case, $\lambda_{0}$. Therefore, the bending length can be approximated as $\ell_{\rm b}\approx r\approx(L/c)-\lambda_{\rm max}\approx A(1-\alpha\varepsilon)$, where $A$ and $\alpha$ are constants determined by the material properties of the yarn. The bending energy scales as $E_{\rm bend}\sim\ell_{\rm b}^{-1}\sim(1-\alpha\varepsilon)^{-1}$. This implies that the high-stress regime scales as $\sigma_{\rm high}\sim\partial E_{\rm bend}/\partial\varepsilon\sim(1-\alpha\varepsilon)^{-2}$. This is consistent with our elastica analysis in the preceding section, where the high-tension limit $q\gg 1$ recovers the same bending energy scaling form.

While the low-stress regime is determined by topology, the high-stress regime is dominated by the material properties of the yarn. Combining these limiting behaviors leads us to a stress-strain relationship $\sigma(\varepsilon)\approx\sigma_{\rm low}(\varepsilon)+\sigma_{\rm high}(\varepsilon)$, which fits our experimental and simulation data quite well, as shown in Fig. 2, Supplementary Fig. 16, and Supplementary Fig. 6.

Our constitutive model for a sample of knitted fabric under uniaxial stress is given by

$$\left\{\begin{array}[]{c}\sigma_{xx}(\varepsilon_{xx},\varepsilon_{yy})=C^{0}_{xxxx}\varepsilon_{xx}+C^{0}_{xxyy}\varepsilon_{yy}+\beta_{xx}\left(\frac{1}{(1-\alpha_{xx}\varepsilon_{xx})^{2}}-1-2\alpha_{xx}\varepsilon_{xx}\right)\\[10.0pt] \sigma_{yy}(\varepsilon_{xx},\varepsilon_{yy})=C^{0}_{yyyy}\varepsilon_{yy}+C^{0}_{yyxx}\varepsilon_{xx}+\beta_{yy}\left(\frac{1}{(1-\alpha_{yy}\varepsilon_{yy})^{2}}-1-2\alpha_{yy}\varepsilon_{yy}\right)\end{array}\right.\,,$$

(Supplementary Equation 14)

where $C_{xxxx}^{0}$, $C_{yyyy}^{0}$, $C_{xxyy}^{0}$, and $C_{yyxx}^{0}$ are components of the rigidity tensor $C_{ijkl}\equiv\partial\sigma_{ij}/\partial\varepsilon_{kl}$, evaluated in the low-strain, linear elastic limit wherein $\sigma_{ij}\approx C_{ijkl}^{0}\varepsilon_{kl}$. The constants $\alpha_{xx}$ and $\alpha_{yy}$ characterize the finite extensibility of knitted fabric when stretched in orthogonal directions, with $\beta_{xx}$ and $\beta_{yy}$ setting the stress scale of the strain-stiffening regime. Our experimental results find significant asymmetry between $C_{xxyy}^{0}$ and $C_{yyxx}^{0}$ components, see Supplementary Table 9, Supplementary Table 11, and Supplementary Table 19. This is backed up by simulations showing significant asymmetry between $x$-strain and $y$-strain. Therefore, we do not enforce the standard symmetry $C_{xxyy}^{0}=C_{yyxx}^{0}$ in our model. We hypothesize that this is due to changes in non-local contact interactions that occur when the fabric is strained in different directions.

We obtained values of these eight parameters for each fabric by fitting the constitutive relations to data via a least-squares scheme. For each fabric, we obtained two independent data series from the uniaxial stretching experiments: (i) measured triplets $S_{x}=\left\{\left(\varepsilon_{xx}^{\rm data},\varepsilon_{yy}^{\rm data},\sigma_{xx}^{\rm data}\right)\right\}$ of strain and stress in the x-direction and (ii) measured triplets $S_{y}=\left\{\left(\varepsilon_{yy}^{\rm data},\varepsilon_{xx}^{\rm data},\sigma_{yy}^{\rm data}\right)\right\}$ of strain and stress in the y-direction. We then minimize the functional

$$I=\mkern-36.0mu\sum_{(\varepsilon_{xx}^{\rm data},\varepsilon_{yy}^{\rm data},\sigma_{xx}^{\rm data})\in S_{x}}\mkern-36.0mu\left(\sigma_{xx}(\varepsilon_{xx}^{\rm data},\varepsilon_{yy}^{\rm data})-\sigma_{xx}^{\rm data}\right)^{2}+\mkern-36.0mu\sum_{(\varepsilon_{yy}^{\rm data},\varepsilon_{xx}^{\rm data},\sigma_{yy}^{\rm data})\in S_{y}}\mkern-36.0mu\left(\sigma_{yy}(\varepsilon_{xx}^{\rm data},\varepsilon_{yy}^{\rm data})-\sigma_{yy}^{\rm data}\right)^{2}$$

(Supplementary Equation 15)

with respect to the seven unknown parameters in the constitutive relation. However, minimizing this functional alone is insufficient because it ignores constraints imposed by the boundary conditions of the fabric. The stress-free boundary conditions couple longitudinal and transverse strains, giving rise to the Poisson effect. In order to introduce this coupling when fitting the data, we determine the pair of Poisson ratios, $\nu_{yx}$ and $\nu_{xy}$, via linear fits to data sets $S_{x}$ and $S_{y}$, respectively. Then we minimize the functional $I$ under the constraints that the linear rigidity tensor components are consistent with these Poisson ratios via the relationships $\nu_{yx}=C^{0}_{yyxx}/C^{0}_{yyyy}$ and $\nu_{xy}=C^{0}_{xxyy}/C^{0}_{xxxx}$. Values obtained for the fitting parameters are shown in Supplementary Table 9, Supplementary Table 11, and Supplementary Table 19. The Young’s moduli are then given by

$\displaystyle\begin{split}Y_{x}&=C^{0}_{xxxx}-\frac{C^{0}_{xxyy}C^{0}_{yyxx}}{C^{0}_{yyyy}}=(1-\nu_{xy}\nu_{yx})C^{0}_{xxxx}\\ Y_{y}&=C^{0}_{yyyy}-\frac{C^{0}_{yyxx}C^{0}_{xxyy}}{C^{0}_{xxxx}}=(1-\nu_{yx}\nu_{xy})C^{0}_{yyyy}\end{split}$

(Supplementary Equation 16)

where we see that $Y_{x}/Y_{y}=C^{0}_{xxxx}/C^{0}_{yyyy}$, yet $Y_{x}\neq C^{0}_{xxxx}$ and $Y_{y}\neq C^{0}_{yyyy}$. Values of the Young’s moduli, along with error estimates based on the variance obtained from least squares fitting, are shown in Supplementary Table 10 and Supplementary Table 12. These figures are plotted in Supplementary Fig. 15.