Knitting is one of the most efficient and widely used techniques for producing fabric membranes. The recent advances in computational knitting make it possible to produce large complex three-dimensional surfaces in one piece without seams [26, 37]. On a modern programmable flat-bed knitting machine it is possible to produce even non-developable surfaces by a local variation of the stitch pattern consisting of an interlooped yarn. Most promisingly, the yarn can be replaced or integrated with electroactive or conductive yarns to produce novel interactive textiles with sensing and/or actuation capabilities [24, 25, 34]. Knitted membranes are usually very flexible because the primary deformation mechanism for the yarn is bending rather than axial stretching. Their unique stretchability and drapability properties make knitted membranes appealing as a reinforcement in composite components [21, 14]. If needed, the stiffness can be increased by inserting straight high-strength fibres during the knitting process. Such reinforced membranes have been recently used as a formwork in architectural engineering . Considering the recent advances in knitting, there is a need for efficient computational approaches for the analysis of large-scale knitted membranes.
For large-scale analysis of knitted membranes, the computational homogenisation approaches which take into account the deformation of the interlooped yarn on the microscale are crucial. There is an extensive amount of literature on finite element-based computational homogenisation of heterogeneous solids, see e.g. the reviews [33, 13, 39]. Two-scale homogenisation often referred to as FE, combined with a yarn-level and a membrane-level finite element model has been also applied to woven and knitted fabrics [29, 12, 22, 11]
. The boundary conditions of the microscale representative volume element (RVE) are given by the membrane deformation and in turn the averaged yarn stress in the RVE yields the membrane stress. However, such schemes are inefficient for knitted membranes with large deformations because of the need to solve a nonlinear problem at each quadrature point of the membrane. It is increasingly apparent that for nonlinear problems two-scale homogenisation must be considered in combination with a data-driven machine learning model
. The model can be trained in an off-line learning phase by solving the microscale problem for a sufficiently rich set of deformation states. Subsequently, the trained model provides a closed-form constitutive equation that is used in the macroscale model. As a machine learning model, for instance, neural networks or GPR models have been used[2, 20, 31, 51]. Alternatively, it is possible to formulate the macroscale finite element problem directly on the training data set bypassing the need to train a machine learning model [17, 35]
. GPR regression is a well-studied Bayesian statistical method and provides as such a principled approach to dealing with epistemic and aleatoric uncertainties and issues such as overfitting. Therefore, we approximate the response of the microscale yarn model with GPR. Their limitation to relatively small dimensions and number of training points is not relevant for this paper.
Owing to the relative slenderness of the yarn on the microscale and the membrane on the macroscale, they are best modelled as a rod and a shell, respectively. Evidently, this leads in comparison to a 3D solid model to an immense reduction in the number of degrees of freedom. The geometrically exact rod theories pioneered by Simo et al.[43, 42] provide a consistent and efficient framework for modelling rods undergoing finite deformations. Later contributions to geometrically exact rod theories include [3, 19, 16]. While in most of these classical rod models the transverse shear deformations are taken into account, in more recent Euler-Bernoulli type models they are neglected [27, 1]. The omission of the transverse shear effectively sidesteps the shear locking problem which is a major impediment in the analysis of slender rods. These new models are usually discretised using smooth B-spline basis functions due to the presence of the higher-order displacement derivatives in the energy functional. The Euler-Bernoulli type rod model introduced in this paper takes into account the stretching, bending and torsion of the yarn as well as the non-frictional rod-to-rod contact between yarns. Similar to the yarn, the fabric membrane can be modelled with geometrically exact shell theories going back again to Simo et al. [40, 41]. As for rods, in more recent Kirchhoff-Love type shell models the transverse shear deformations are neglected and the weak form is discretised using smooth B-splines or subdivision basis functions [7, 8, 18]. Although the bending resistance of knitted membranes is usually very small, a suitably chosen shell bending energy term can serve as a regularisation for wrinkling under compressive stresses [8, 5]. In this paper we make use of the Kirchhoff-Love subdivision shell implementation introduced in [6, 23].
The outline of the paper is as follows. In Section 2 we introduce our finite deformation rod model, its discretisation with B-splines as well as the treatment of rod-to-rod contact. Subsequently, in Section 3 we discuss the proposed microscale yarn-level RVE model for computational homogenisation. This model is verified and validated with experimental and numerical results from literature. In Section 4, we introduce the data-driven GPR model and describe its training with the yarn-level RVE model. Finally, in Section 5, we first discuss the training of the GPR model and then analyse membranes subjected to tension to assess the accuracy of the obtained data-driven constitutive model.
2 Finite deformation analysis of yarns
In this section we summarise the governing equations for the finite deformation Euler-Bernoulli rod model for the yarn. The presented equations are without loss of generality restricted to rods with circular cross-sections. We take into account rod-to-rod contact by enforcing the non-penetration constant with the Lagrange multiplier method.
The geometry of the rod is described by a set of circular cross-sections connected by their line of centroids. In accordance with the Euler-Bernoulli assumption, equivalent to the Kirchhoff-Love assumption for shells and plates, transverse shear is neglected so that cross-sections remain always normal to the line of centroids.
The position vectors of material points in the reference and deformed configurations and are parametrised in terms of the convective coordinates as
where and denote the lines of centroids parameterised by , see Figure 1. In turn, the cross-section with the radius is parameterised by and . The tangent vectors to the line of centroids are given by
The two orthonormal directors and are chosen so that they satisfy in the reference configuration
Here and in the following Greek indices take the values and summation over repeated indices is assumed. The two orthonormal directors in the deformed configuration are obtained by rotating the reference configuration directors with a rotation matrix ,
This, in combination with the Euler-Bernoulli assumption, ensures that the two directors and in the deformed configuration satisfy
The rotation is composed of two rotations,
That is, the reference directors are mapped to the deformed directors in two steps using an intermediary configuration with . The matrix describes a rotation by an angle about the unit tangent and is according to the Rodrigues formula given by
with the identity matrix
and the skew-symmetric matrix
To derive the strains corresponding to the assumed kinematics (1
), we consider the Green-Lagrange strain tensor of 3D elasticity
Here and in the following Latin indices take the values . The covariant basis vectors and and the contravariant basis vectors and are defined as
where is the Kronecker delta. The corresponding two metric tensors and are given by
We identify as the membrane, bending about , bending about and torsional shear strains, respectively. Moreover, due to the Euler-Bernoulli assumption the strain components and are zero and the in-plane shear strains and are induced only by torsional shear strain .
2.2 Equilibrium equations in weak form
The potential energy of a rod with the line of centroids and the cross-section occupying the volume in its reference configuration takes the form
where is the strain energy density and is the potential of the externally applied forces. At equilibrium the potential energy of the rod is stationary, i.e.,
|with the external virtual work and the internal virtual work|
where is the second Piola-Kirchhoff stress tensor. The strain tensor of the rod (13) depends on the displacement of the line of centroids
and the rotation angle . Hence, we can write for the internal virtual work (15b) more succinctly
where and are the internal forces conjugate to the virtual displacements and rotations .
As a material model we use the isotropic St Venant-Kirchhoff model with the strain energy density
and the fourth-order constitutive tensor
where and are the two Lamé parameters . The contravariant metric tensor is determined from the relation .
To derive analytical expressions for the internal forces we first introduce the rod strain tensor (13) and the constitutive equation (19) in the internal virtual work (17). Subsequently, we integrate over the rod cross-section analytically to obtain the internal forces
Here, the axial force
, the bending momentsand the torque are defined as
where , and are the cross-section area, second moment of area and torsional constant of the circular rod, and is its Young’s modulus. The tedious but straightforward derivation of internal forces and their derivatives are summarised in Appendix A.
2.3 Finite element discretisation
We follow the isogeometric analysis paradigm and use univariate cubic B-splines to discretise the lines of centroids and in the reference and deformed configurations. We choose smooth B-splines because the bending strains in (13) require at least continuous smooth basis functions. The kinematic relationship (16) is restated after discretisation as
where the B-spline basis and their coefficients correspond to the control vertices on the discretised rod centreline.
2.4 Yarn-to-yarn contact
We use the Lagrange multiplier method to consider pointwise non-frictional contact between two circular rods. By closely following Wriggers et al.  and Weeger et al. , we add to the total potential energy in (14) the contact potential energy
where the Lagrange multiplier represents the repulsive normal force between the two rods, and the non-positive gap function depends on the minimum distance between the two rods. The distance between the lines of centroids of the two rods and is given by
and its minimum by
This minimum can be determined using the Newton-Raphson scheme. Ultimately, the non-positive gap function is defined as
where is the radius of the rods.
2.5 Verification of the rod model
We consider a membrane-bending-torsion interaction problem to verify the accuracy of the presented rod formulation. A spatial helicoidal spring is clamped at one end and a vertical load of magnitude 25 kN is applied at the other end, see Figure (a)a. The reference geometry of the line of the centroids is given by
The material and geometric properties of the spring including the reference director orientations are given in Figure (a)a. The tip load is increased from zero to 25 kN in 10 uniform load steps. In Figure (b)b the deflected shapes for five different load levels are shown. As can be seen in Figure 6 the obtained tip displacements are in excellent agreement with the results presented in Bauer et al. , in which a slightly different approach was used to parameterise the rotations of the two rod directors .
3 Computational Homogenisation
3.1 Microscale analysis
A characteristic representative volume element (RVE) as depicted in Figure 7 is chosen to represent the periodic microstructure of a weft-knitted membrane. The intricate spatial arrangement of the yarns in the RVE is modelled, as in recent works [10, 49], using the approximate geometry proposed by Vassiliadis et al. , see Appendix B. For alternative geometric descriptions see  and references therein. First the yarn centrelines are defined and then the yarns are generated by sweeping a circle along those lines.
On the macroscale we model the membrane with subdivision shell finite elements [8, 6, 23] and consider for homogenisation only the in-plane membrane response. The very small out-of-plane bending stiffness contribution of the rods to the membrane bending stiffness is not taken into account.
In first-order homogenisation, the macroscale membrane deformation gradient is used to define the boundary conditions on the RVE. Here and in the following the macroscale and microscale quantities are denoted by subscripts and , respectively. It is assumed that the volume averaged deformation gradient of the rod in the RVE and the membrane deformation gradient are equal, that is,
where is the RVE volume and is the rod volume within the RVE in the reference configuration. Without going into details, the deformation gradient of the line of centroids is given by the kinematic assumption (1). The Hill-Mandel lemma states that the macroscale and microscale work for an RVE must be equal
where is the macroscale first Piola-Kirchhoff membrane stress and is the microscale first Piola-Kirchhoff rod stress. The integral represents the internal virtual work of the rod within the RVE and is given by (17). For the lemma (30) to hold only certain types of RVE boundary conditions can be chosen, including Dirichlet and periodic; see [13, 52] for details. Moreover, as a consequence of this lemma the macroscopic first Piola-Kirchhoff stress can be obtained from the volume averaged internal energy density
where the set denotes the deformation gradients satisfying the RVE boundary conditions. The microscopic energy density of a rod with an isotropic St Venant-Kirchhoff material is given by (18).
In this work we choose as RVE boundary conditions the periodic boundary conditions depicted in Figure 10. Due to the orthotropy of the RVE response, the biaxial stretching and shear response are decoupled. Therefore, we consider the two cases separately and choose for each different boundary conditions. Boundary displacements in the thickness direction are zero to simulate plane stress conditions.
3.2 Verification and validation of the microscale model
To verify and validate our microscale model we consider RVEs under wale-wise and coarse-wise uniaxial tension and shear as reported in Dinh et al.  and Weeger et al. . For a detailed problem description we refer to  and . In our computations the yarn is discretised with 128 rod finite elements. In Dinh et al.  the yarn is modelled as a 3D solid and in Weeger et al.  as a discretised Euler-Bernoulli rod with the collocation method. As evident from Figure 11 our results are in excellent agreement with the experimental and computational results obtained with other approaches. Throughout this paper, we use the yarn geometry and material parameters given in the above-mentioned two papers.
4 Data-driven homogenisation
4.1 Review of Gaussian process regression
Gaussian process regression is a statistical inference method rooted in Bayesian statistical learning. In the following, we briefly review the key steps in Gaussian process regression. For further details we refer to Rasmussen et al. . To begin with, we assume as a prior in the Bayesian formulation a random function given by the zero-mean Gaussian process
The respective covariance function is chosen to be the (stationary and isotropic) squared exponential function
are the yet to be determined scaling and characteristic lengthscale hyperparameters. For convenience, we defineas the vector of hyperparameters. The infinitely smooth covariance function (33) encodes our prior assumptions about before observing any training data. Other choices are possible, see [38, Chapter 4].
Next, we consider a database comprised of a known training dataset and testing points. Each data point is given by an input vector and a corresponding scalar observation , with . All training data points are collected in (, ) and all testing data points with unknown in (, ). Moreover, we define the covariance matrix with the components . The discretisation of the Gaussian process (32
) for the considered training and test data is given by the multivariate Gaussian distribution
The probability density of the unknown dataat the given locations conditioned on the known training data reads
Hence, the best estimate foris given by the expectation
and the uncertainty of the estimate is given by the covariance matrix
The covariance matrix of size is dense for the covariance kernel (33). If the number of training points become too large to invert , alternative approaches leading to sparse covariance matrices need to be considered [38, Chapter 8].
To estimate the hyperparameters , we consider the marginal likelihood, or the evidence,
with the likelihood and the prior , where are function values evaluated at , c.f. (32). The marginal likelihood is the probability of observing the data for given hyperparameters . Hence, it suggests itself to choose the hyperparameters so that the probability observing the given training data is maximised. In practice, it is numerically more stable to compute the maximum of the log marginal likelihood given by
Note that the is a monotonic function so that
We use the scikit-learn Python library  to find with a gradient-based optimisation algorithm. The marginal likelihood is a non-convex function and certain care has to be taken to find the global maximum. After the optimal hyperparameters are determined they are used in computing the mean (36) and covariance (37).
4.2 Gaussian process homogenisation
The implemented data-driven homogenisation framework is illustrated in Figure 12 and closely follows Bessa et al. . Data-driven homogenisation begins with constructing a response database for the RVE. This response database is designed by defining the microscale boundary value problem (BVP) using three types of design variables, namely geometric properties , material properties and boundary conditions . In this paper, we concentrate on the design variable relevant to homogenisation. It is straightforward to include and which lead to the material design of knitted membranes.
Next, we determine the hyperparameters of the GPR model using -fold cross-validation rather than directly using the obtained with (40) considering the entire training dataset. The GPR model is iteratively evaluated by using folds for training and the remaining one for testing. We choose which provides a good compromise between accuracy and efficiency. As error metrics for the predicted strain energy and stress resultants we use the correlation of determination and the mean squared error (MSE) given by
where is the number of the test points,
is the target output,is the predicted output and is the expectation of the target output given by the GPR model. When evaluating the MSE of the predicted stress tensor, we take the square root of the squared sum of the MSE of each stress resultant
The GPR model training is implemented in a Python environment using scikit-learn machine learning library . For each test fold, we determine the corresponding hyperparameters by maximising the log marginal likelihood (40). It is expected that the obtained hyperparameters have similar values for all test folds. In practice, this is not the case because the log marginal likelihood is usually a non-convex function. Therefore, we vary during optimisation the upper and lower limits of the hyperparameters so that all test folds yield similar hyperparameter values. Once this is achieved, we record the hyperparameters of the GPR model, which has the lowest MSE, and save the model for subsequent macroscale analyses.
Lastly, we integrate our in-house thin-shell solver [7, 23] with the trained GPR model to simulate the homogenised knitted membrane. As shown in Figure 12, for given macroscale deformations (strains) the trained GPR model yields the corresponding predicted stress resultants and their tangents to the macroscale shell solver.
4.3 Macroscale analysis
For plane-stress membrane deformation, a response database is created with the three in-plane components of the Green-Lagrange strain as the inputs and RVE volume averaged strain energy as the output, i.e. and in Section 4.2. With a slight abuse of notation is here, in contrast to Section 2.1, a two-dimensional second order tensor. In the four-dimensional response database, points have the coordinates . The strain components and correspond to a biaxial deformation state that is decoupled from the shear deformation state with , see Section 3.1. The RVE strain energy is obtained by considering a biaxial strain state with and and a shear deformation state with and summing up their strain energies. We use the validated strain limits of Figure 11 in constructing the response database. Thus, allowing for a 5 compressive strains, we define the design variables of the response database as,
We solve the RVE problem with the boundary conditions stipulated by the strain states defined in (43) and store the strain components and respective volume averaged energies in the response database. After passing the GPR model training and testing phase, as shown in Figure 12, the strain energy density of the plane-stress deformation is predicted for a given deformed state of the RVE. Plane stresses and stress tangents are computed by differentiating (36), that is,
5.1 GPR model training and testing
We start with presenting the error metrics of GPR model training and testing for varied sizes of response databases. GPR model training errors are observed to take values and . Hence, we present only the error metrics of testing datasets. All the following computations are performed on the Intel Core i5-4590 CPU @ 3.30GHz 4 processor.
Figure 13 shows the MSEs of predicted energy density and second Piola-Kirchhoff stress resultants, as well as a comparison between two sampling techniques, namely uniform sampling and Sobol sampling , on the testing errors. Uniformly sampled data points are generated in a way that for the three features in the input vector , entries occupy the response database, where
is the number of uniformly distributed entries in each feature. In total seven datasets withare considered; of each are taken as training points and the remaining as testing points. Moreover, is rounded to the fifth decimal number for all response databases. MSEs in stress predictions are comparatively higher than those in energy predictions but remain less than 2. Considering the error convergence and algorithmic efficiency, we use a response database with 2601 test data points in the subsequent macroscale simulations. Furthermore, this particular trained GPR model has the optimum hyperparameter values and the maximum log marginal likelihood . Model training time was recorded as 29 minutes and 48 seconds.
In Figure 16, we visualise and compare the predictions of a subsample of the chosen response database. Figure (a)a depicts the strain energy density predictions, whereas Figure (b)b presents the stress predictions for the biaxial homogenised response of an RVE.
5.2 Stretched membrane I: comparison of yarn-level and homogenised displacements
We consider a membrane under uniaxial tension and compare its response using the Gaussian process homogenised material model and an equivalent yarn-level model. The yarn-level model comprises 12 loops in the course direction and 12 loops in the wale direction. We use for the analysis with the homogenised material model a membrane of equal size, that is, long, wide and thick. The membrane is discretised with a structured quadrilateral mesh with elements. Problem descriptions of the two models are presented in Figure 19.
Both models are stretched by in the -direction and the resulting and displacements are presented in Figure 22. For comparison purposes, the results of the yarn level model in Figure (a)a are visualised by projecting the nodal values of the yarn onto the
plane and interpolating them on a Delaunay mesh. The-direction displacement distribution of the yarn model in Figure (a)a is slightly different from that of the homogenised membrane in Figure (b)b. This difference has two causes. First, zero out-of-plane displacement boundary conditions on the top and bottom edges of the yarn-level model directly contribute to the observed difference. These boundary conditions are applied to simulate selvedge stitches that restrict the top and bottom yarns from being freely straightened during deformation. Secondly, the geometric asymmetry of the yarn-level model about the mid-horizontal plane has an influence on the overall response due to the relatively small number of loops used in both course and wale directions. Furthermore, a very similar Poisson effect is observed in both the yarn-level and homogenised membrane models. The absolute maximum -displacement is recorded with for a stretch by in the -direction.
5.3 Stretched membrane II: homogenised stresses
We perform two stretch tests on a square knitted membrane sheet with a side length of in the course and wale directions. A structured quadrilateral mesh of size ( elements) is used. During the displacement controlled deformation one edge is fixed while the other is stretched in the course or wale direction, respectively, by . The deformed shapes and stress contours are depicted in Figure (a)a for course-wise and in Figure (b)b for wale-wise stretching. Figure 25 clearly shows the orthotropic response of the knitted membrane as the stress resultants are different depending on the direction of the stretching, comparing, e.g., in Figure (a)a with in Figure (b)b. Furthermore, the stiffer response in the wale-wise direction manifests itself in a higher stress in Figure (b)b in comparison to in Figure (a)a.
We introduced a data-driven approach for computational homogenisation of knitted membranes to address the challenges posed by conventional homogenisation schemes. The inherent large deformations in knitted textiles are accurately captured by the finite deformation rod model and validated against experimental and numerical results. Incorporating the statistical GPR model in computational homogenisation circumvents the need for expensive microscale simulations at each quadrature point of the macroscale membrane, thus yielding significant computational savings.
The presented approach can be extended in several ways. In this paper, we considered only weft-knitted membranes with a uniform stitch pattern. To produce complex three-dimensional surfaces, it is necessary to alter the stitch pattern, for instance, by locally increasing or decreasing the number of loops from row to row or reducing the number of rows. Such changes represent discontinuities in the stitch pattern, and suitable RVEs have to be defined to capture their homogenised response. Putting aside questions of RVE size and validity of homogenisation assumptions, it is straightforward to consider in the introduced data-driven approach RVEs with other stitch patterns. Furthermore, depending on the loading conditions knitted membranes are prone to geometric instabilities in form of wrinkling on the macroscale and rod buckling on the microscale. In presence of such instabilities the homogenisation assumptions are usually not valid and alternative approaches must be used [30, 28]. The data-driven variants of these approaches are essential for the analysis of large-scale knitted membranes. Finally, a key advantage of the data-driven approach is the prospect to consider design parameters, pertaining to the stitch geometry or yarn material, in the GPR model. This opens up the possibility to optimise those parameters with an efficient gradient-based optimisation algorithm.
Appendix A Derivatives of strains and internal forces
In this Appendix we summarise the detailed equations used in the implementation of the introduced finite deformation rod finite element.
a.1 Strain derivatives
The derivative of the components of the strain tensor (13) with respect to the nodal displacements are given by
The derivative of the rotation matrix in (9) with respect to the nodal displacements reads
The derivatives of the bending and torsional shear strains in (13) with respect to the nodal twist are given by