Mechanical modelling of 3D woven composites considering realistic unit cell geometry.

Modelling the mechanical performance of textile composites is typically based on idealised unit cell geometry. However, 3D woven composites feature more complex textile architecture then 2D woven materials, and in reality nominally straight warp and weft yarns can also possess signiﬁcant waviness. For such textiles, idealising yarns as straight entities becomes an oversimpliﬁcation. In this study, the voxel method and a continuum damage model are used in a ﬁnite element analysis to compute stress–strain curves for an orthogonal 3D woven composite under tensile loading. The main goal of this study was to compare results produced using idealised geometry with realistic geometry obtained from detailed simulation of the preform during weaving and compaction. Signiﬁcant variation in predictions was obtained using the different geometrical models. The idealised model lead to an overestimation of stiffness and strength compared to experiment due to the neglecting of yarn waviness, whereas the simulated geometry models produced more conservative results closer to experiment. 2014 The Authors. Published


Introduction
3D woven composites can provide a potential solution to the fundamental limitations of traditional laminated composites; delamination and labour intensive manufacture. The addition of binder yarns provides through-thickness reinforcement leading to greatly enhanced interlaminar properties as well as binding the fabric to enable near-net-shape preforms to be woven and handled. However, despite these advantages, the use of 3D woven composites has been largely limited to niche applications. One of the key reasons for this is the lack of predictive numerical tools, which limits their ability to be used at the early stages of design.
Mechanical performance modelling of textile composites typically begins with the definition of textile unit cell geometry using a specialist pre-processor such as TexGen [1,2] or WiseTex [3,4]. Such software can produce an idealised representation of a textile, sufficient for the modelling of many types of composites with 2D reinforcement. However, some 3D woven textiles present a significant challenge to model due to their inherent complexity. While TexGen is capable of modelling such textiles, the idealised geometries that it produces can neglect realistic features such as yarn waviness and yarn pinching which play a significant role in determining their resulting properties, especially strength [5,6]. Previous work [7] used a finite element (FE) model, based on the digital element method [8] in order to predict the deformations of an orthogonal 3D woven fabric at the unit cell level. Each yarn was represented as a bundle of 61 chains of beam elements, with contact between chains within a yarn defining yarn cross-sectional shape and contact between each yarn assembly defining the yarn paths. Tensioning of the binder yarns was used to simulate the weaving process prior to moulding through compaction between rigid plates. Throughout this procedure the unit cell was subjected to periodic boundary conditions such that the resulting geometry was periodic.
Conventional FE models of 2D woven textile composites are typically meshed with tetrahedral elements, with the matrix domain designated as the inverse of the yarn volume within the global unit cell domain (see for example [9]). Nodes of matrix and yarn elements at the interface between the two different material phases must be shared to form a conformal mesh. However, realistic models of complex 3D woven architectures will likely feature degenerated resin regions and small yarn interpenetrations meaning that such meshing practices are no longer feasible, especially when a periodic mesh is required. Other methods have been proposed to model such textile composites including the domain superposition technique (DST) [10] and embedded element technique [11]. In these methods, the matrix domain is assigned as the global domain leaving the yarns to be meshed independently. A  ensures continuity of displacements while properties of the yarns are reduced to account for the superposition. The independent mesh method (IMM) [12] is a similar technique where any matrix elements lying fully within yarn domains are excluded from the analysis. Matrix elements which lie on the boundary are then subjected to a refinement scheme. The extended finite element method (XFEM), originally developed for modelling discontinuities due to cracks, has also been applied to heterogeneous materials with discontinuities at the material interface [13], including textile composites [14]. The global domain is meshed using regular hexahedral elements and the nodes of elements in the interface region are enriched with extra degrees of freedom to account for the strain field jump. However, it is still required that the geometry is free from interpenetrations. In this study, a technique based on voxels [15]; meaning a 3D pixel, is used. The global domain is meshed in a similar process to XFEM and elements are assigned to either the matrix or yarn domain depending on the location of the element centroid. The voxel method has been shown to be adequate for mechanical modelling of 2D textile composites [16].

Textile geometry
The accuracy of mechanical modelling of textile composites for strength prediction is dependent on the modelling technique, material models and assumed textile geometry. A previous project considering failure modelling of 3D woven composites [17] concluded that a significant limitation was with regard to the idealised geometries being used. For textiles with complex internal architecture, there was no further merit in pursuing relatively small improvements with the failure models whilst significant discrepancies existed between considered and real geometries. A schematic representation of an orthogonal 3D woven fabric is shown in Fig. 1. The fabric has two sets of binder yarns, each arranged in a 5 harness satin style. One set of these binder yarns float on the upper surface of the fabric while the other float on the bottom surface. The original idealised TexGen unit cell model of this fabric used in the aforementioned study is shown in Fig. 2(a). In this geometrical model, the in-plane warp and weft yarns are completely straight with zero crimp, waviness or cross-sectional variation. Binder yarns follow a path with straight horizontal or vertical sections, constant cross-section and extend beyond the surface of the fabric, resulting in a resin layer on the surfaces of the composite. None of these features are representative in the real moulded composite.
However, there have been recent developments to TexGen aimed specifically at improving the quality of 3D woven textile geometry. Despite this, it was not possible to generate a model at the level of compaction consistent with the real infused composite (58.5% overall fibre volume fraction (o-FVF)) using the latest version of TexGen (3.5.3). Fig. 3 shows an attempt at generating such a model. Fig. 3(a) shows excessive squeezing of binder yarns at the surface while Fig. 3(b) shows penetrations of binder and weft yarns. A valid model at the target o-FVF could only be achieved by shrinking yarn cross-sections, thus requiring excessive intra-yarn fibre volume fraction (iy-FVF). As such, the new idealised geometrical model used in this study has a lower o-FVF of 52.0%, similar to the original idealised model, and is shown in Fig. 2(b). This model features binder yarns which are flush with the fabric surface, achieved using binder yarns which are near circular in the vertical sections but flatten considerably on surface while inducing crimp in the surface weft yarns. However, elsewhere the in-plane yarns remain straight with constant cross-section.
A numerical model has been developed and used to predict the deformations of the 3D woven fabric considered in this study during weaving and subsequent compaction in a mould tool [7]. This was conducted in order to produce accurate geometry for unit cell analysis of textile composites. Whilst these deformation models can produce geometry at levels of compaction as high as the real fabric, for comparative purposes, this simulated geometry TexGen model is also shown at 52.0% o-FVF in Fig. 2(c), and is visually considerably different to the idealised models. It features significant in-plane and out-of-plane waviness in warp and weft yarns and variation of cross-sections in all yarns.
The average waviness levels in each of the geometrical models are listed in Table 1. This was determined by splitting the yarn path into short segments and calculating the angle of each segment from the nominal path both in-plane (h IP ) and out-of-plane (h OOP ), as shown in Fig. 4 and Eq. (1). The overall waviness (h overall ) of each segment was assessed by considering a 3D vector, a, defined by h IP and h OOP . h overall was calculated as the angle of a with another vector, b, defining nominal yarn path, as calculated in Eq. (2). The final waviness value is the average of this value over the total length of the yarn.
where : The output of the deformation model has been validated against computed tomography (CT) scans of the dry fabric subjected to insitu compaction across a range of volume fractions in [7], with excellent qualitative correlation. A CT image of a unit cell is presented in Fig. 2(d) where it can clearly be seen that the simulated geometry model represents the architecture of the real fabric much more accurately than the idealised models, but unfortunately experimental data for the waviness of the fabric at the high level of compaction considered in this paper is not available.

Conversion of realistic geometry from deformation models
The output of the deformation models, with yarns represented as a bundle of beam elements, is not immediately useful for meso-scale mechanical performance analysis where yarns require definition as solid geometry. Conversion of the deformation model geometry was performed using Python scripting to create the simulated geometry TexGen model [18] shown in Fig. 2(c), from which FE models can be generated. Building a unit cell textile model in TexGen requires definition of three attributes; yarn paths, yarn cross-sections as well as yarn repeats with a domain. For detailed information regarding the theory of geometrical modelling of textiles as applied in TexGen, refer to [1].
The path of a yarn in TexGen is represented by a one-dimensional line defined in three-dimensional space. The script was used to march along each yarn, stopping at regularly spaced intervals to  calculate the yarn centroid and attribute these as master nodes which were joined with a periodic spline (Fig. 5). Yarn crosssections are defined in TexGen as two-dimensional shapes in a plane perpendicular to the yarn tangent. Planes were specified at each master node to account for the variation in cross-sectional shape along the yarn path and the intersection points of each beam element within the yarn was mapped to these planes. Next, cross-sections were defined based on these mapped points, using the procedure outlined below to generate polygons. These polygons are sampled by a number of equally spaced slave nodes for graphical rendering and meshing. Convex hull algorithms are widely used in computational geometry to define the smallest convex polygon around a set of points e.g. [19,20]. This can be visualised using an analogy in which the points are represented by pins in a board. By tying a piece of string to a starting pin, the convex hull is the shape formed by wrapping the string around the set of the pins until it returns to the starting pin. However, the shape of a yarn cross-section is not always convex and applying the convex hull would give an inaccurate representation of the yarn in such cases. Therefore a modified convex hull algorithm was implemented, beginning by defining the first point in the polygon as the node with the largest x co-ordinate (in an x-y plane), then proceeding anti-clockwise around the section. Using the previous analogy, a string of finite length, r, is employed to define a search radius for the subsequent selection which is the point creating a line with the smallest external angle from the previous polygon section. A closed polygon is formed once the initial node is reselected and is subsequently expanded by the radius of the beams, d/2. The search radius, r, was defined as r = xd, where the parameter x can be used to adjust the level of smoothing of the cross-section. A value of x = 4 was found to produce good results and is compared with the convex hull solution for two extreme examples of yarn cross sections generated by the deformation models in Fig. 6. While the convex hull produces good results for the thick yarn, its representation of the thinner yarn is poor since the large concave regions are smoothed over. The modified algorithm, however, produces good results in both situations.
Since the output of the unit cell deformation models is periodic, specification of yarn repeat vectors corresponding to the tessellation of the unit cell can be utilised in TexGen to define an infinite fabric which is then trimmed to a finite domain of the unit cell.

Modelling methodology
In this study, tensile loading in warp and weft directions will be considered for the following four scenarios: A. New idealised geometry model (52.0% o-FVF) B. Simulated geometry model (52.0% o-FVF) C. Simulated geometry model (58.5% o-FVF) D. Experimental data (58.5% o-FVF) The quoted o-FVFs for A-C refer to the TexGen geometrical models. The voxel models generated from these TexGen models show negligible difference from these values, since any errors in yarn volume fraction (YVF) resulting from the voxel discretisation are offset by adjustment of the iy-FVF. The results of A and B are used to assess the effect of the different geometries on mechanical performance prediction, since the two models can be directly compared at the same o-FVF. Further insight can be gained regarding the effect of compaction from B to C which can only be achieved by using models of the realistic geometry. The accuracy of each model as compared to experiment will also be assessed.

Mesh and boundary conditions
Voxel meshes can be automatically generated in TexGen with each voxel being designated to either a yarn or the matrix depending on the location of its centroid. This means the technique is tolerant of small yarn interpenetrations in the geometrical models which can be introduced during the conversion process from deformation models. Unit cell models of size 27.8 Â 9.9 Â 5.3 mm were used with cubic shaped voxel elements providing an element length of 0.12 mm. Using this technique, it is trivial to generate periodic meshes of complex geometries thus reducing pre-processing effort and facilitating automation. It also avoids the need to specify the fictitious resin gaps between yarns often resorted to in conventional modelling of textile composites to produce adequate quality elements (see for example [21]). Material orientations were assigned to each voxel element using the procedure applied in TexGen which is outlined in [1]. Periodic boundary conditions for the unit cell with stagger in tessellation shown in Fig. 1(a) were applied as derived in [22]. Since the unit cell represents the full Master node Slave node

Material properties
In each of the models the yarn cross-sectional area and hence iy-FVF varied in segments along yarn path. This variation was from 47-78% iy-FVF in the idealised model, though most of this variation occurs in the binder yarns. The simulated geometry models exhibited a high diversity of yarn shapes and dimensions in all yarns with iy-FVF ranging from 36-87%, where high iy-FVFs occurred in regions of yarn to yarn contact at cross-over points where load is transmitted through fabric thickness during weaving and compaction. Whilst the iy-FVF remains below the maximum hexagonal array fibre packing limit of $91%, since this cannot be violated in such multi filament based deformation models, the high iy-FVFs suggested in certain localised regions is exaggerated. In practice it is difficult to exceed iy-FVFs greater than 70-75% (see for example [23][24][25]) as unlike the model, fibres in a real yarn are not perfectly aligned and have some level of entanglement.
The resin pockets were assigned the properties of the pure isotropic epoxy resin while the yarns were treated as transversely isotropic homogenous material similar to a unidirectional lamina. Properties were applied to each element within the yarn depending on the local iy-FVF, in order to accurately represent the distribution of properties along the yarns. Fibre and matrix properties are shown alongside estimated yarn properties at a nominal iy-FVF of 70% in Table 2. Properties of the carbon fibres and epoxy resin were taken from manufacturer's data, though some fibre properties were taken from the literature, with averaged experimentally derived values used where available.
Yarn elastic constants were calculated from constituents using a micromechanical FE model with hexagonal fibre packing. Strength properties were estimated using a simple analytical model, called the generalised chess model, first used in [26]. Fig. 7(a) shows the micro scale unit cell of a unidirectional fibrous composite with regular square packing. This unit cell has two planes of symmetry and thus the considered region can be reduced to a quarter of this unit cell. Fig. 7(b) shows the quarter unit cell approximated by a square fibre grid consisting of four blocks. Blocks 1, 2 and 4 represent matrix material with block 3 representing the fibre. A set of conditions including the prescribed stress/strain field, continuity of traction on the inter-block boundaries and symmetry form a system of equations relative to the stress in each block. In essence, the obtained representation is similar to the discretisation of the unit cell by four finite elements.
These equations along with the Hooke's law are sufficient to identify the stresses in each block, which are assumed to be uniform. Once the stresses are determined, strength can be estimated. The fibre dominated longitudinal strength is calculated by maximum stress criterion and resin dominated strength properties are calculated using the following criteria: max p mx p cr ; q mx q cr P 1 ð3Þ where P mx and P cr are the maximum hydrostatic tension in the matrix (crucial for failure of brittle epoxies [27,28]) and its critical value. In an elementary uniaxial tensile test the epoxy matrices fail in a brittle fashion and hence p ¼ p cr ¼ S t 3 at failure (one third of the apparent uniaxial strength S t ). q mx and q cr are the Mises stress and the correspondent critical value characterising the probability of a ductile failure. In pure shear the epoxy yields as an ideal plastic material before failing, hence q cr can be taken as the yield stress.
The model was used to predict strength under three uniaxial conditions; transverse tension, in-plane shear, out-of-plane shear. The prediction of the transverse tensile strength follows the anticipated trend; strength drops when iy-FVF is increased from 0% to 40% and remains nearly constant as it further increases. The transverse tension and shear strength values are in good agreement with the available data, e.g. [29].

Material damage models
Carbon fibres often exhibit non-linear tensile behaviour, e.g. at an applied strain of 1% the Young's modulus may increase by 20-25% of its initial value [30]. Recent studies showed this phenomenon (known as the 'non-Hookean' effect) to be partly responsible for the non-linear behaviour of 3D woven composites [31]. However, due to the lack of data for the particular carbon fibres used here, the behaviour of yarns was assumed to be linear prior to damage initiation. The matrix was treated similarly. After damage initiation a continuum damage mechanics (CDM) model [32] was applied both to yarns and matrix. It assumes that damage in a material can be considered as a reduction of stiffness properties over a region of material, a single finite element in this case.
Matrix damage was assessed by the pressure dependent modified von Mises criterion listed in Eq. (4) [33]. Damage initiation was assumed when the damage parameter D m exceeded the value of 1.0.  where r i with i = 1, 2, 3 are the principal stresses and S j m with j = t, c are the matrix tensile and compressive strengths.
Three modes of failure were considered for the yarns: longitudinal tension, transverse tension and shear. The damage state of the material was defined by a set of three parameters D i listed in Eqs. (5)-(7) with indices i = 1, 2, 3 representing longitudinal, shear and transverse failure modes respectively.
where; r ij are components of stress tensor, S i j k are strengths of yarn with indices i, j = 1, 2 corresponding to the directions while k = t, c stands for tension and compression. The moduli of the damaged material were described by Eqs. (8)-(10) for the yarns and Eq. (11) for the matrix: where stiffness properties with a superscript 0 refer to the undamaged value and P(D i ) is a penalty function defined as: It can be seen that damage in the longitudinal direction causes abrupt reduction of longitudinal properties whereas other stiffness properties are subject to an exponential degradation. Poisson's ratios m 12 , m 23 and m 13 remain intact, while the corresponding Poisson's ratios in other directions (i.e. m 21 , m 32 and m 31 ) were adjusted to keep the stiffness matrix of an element symmetric. Parameters c 1 and c 2 in the penalty function are empirical constants where their ratio determines how fast damage will propagate through an element and defines the final failure of the element (Fig. 8). This will have some mesh dependence, with coarse meshes requiring a higher ratio to prevent the dissipated energy increasing. The ratio c 2 /c 1 was initially chosen from a model of a 2D woven unit cell [16] with a more refined mesh. Based on the initial results of loading in the warp direction an improved value of 4.0 was selected. The influence of mesh size whilst keeping this parameter constant has been investigated and results are discussed in Section 4.

Results and discussion
Voxel meshes used for the idealised and simulated geometry models are shown in Fig. 9. It was relatively simple to represent the straight yarns of the idealised model with voxels, but the complex geometry of the simulated geometry models proved more challenging with angled or curved yarns having stepped or jagged surfaces. The quality of the voxel representation of geometry improves with refinement, however the large unit cell size of this fabric placed a practical limit of the level of refinement possible. The baseline models had nearly 800,000 elements and took 5 days to run on a high performance computing system using 32 CPUs.
Young's moduli and ultimate tensile strength values are presented for each of the models compared to experiment in Table 3. Stress-strain curves for each data set are shown for the warp and weft directions in Figs. 10 and 11 respectively. Experimental curves [34] shown are for a representative test and model results using the low compaction geometry have been normalised by o-FVF to aid comparison.
Non-linear behaviour of a 3D woven composite was shown in [31] to be the result of a combination of three phenomena: fibre non-linearity, inelastic yarn straightening and damage development in the composite. The latter of these two phenomena are accounted for in the modelling here. It was also shown in [30] that the first two phenomena dominate at lower strains while damage development overcomes them at higher strains.
The experimental curves in Figs. 10 and 11 are bi-linear with a stiffness reduction kink at around 0.6% strain. It is believed that this non-linearity is predominantly due to damage development,  more specifically transverse yarn, and also resin, cracks. These have been observed in an orthogonal 3D composite of similar weave style but fewer layers in interrupted tests (Fig. 12). Significant damage is also present in the models prior to final failure, see Fig. 13 showing 'cracks' of damaged material forming in the yarns perpendicular to the applied load.
Longitudinal strains in the yarns were assessed in one of the models. It was found that at an applied strain of 0.6% (where $10% of transverse yarns material has sustained some damage) only a small fraction (<1%) of the yarns have strain values higher than 1%. Based on [31], this implies that overall increase of stiffness cannot be higher than 10-15%. It was therefore considered reasonable to neglect fibre non-linearity in the material model.
The result of the idealised model was an overestimation of stiffness in both directions by around 15%. The low compaction simulated geometry model predicted moduli very close to experiment with the high compaction model having slightly reduced values. The same trend is replicated with strength prediction for each model. The damage model used for the analysis is based on a simplified mechanical approach and thus cannot capture transverse failure with great accuracy. Stress distribution, damage onset and final strength are significantly influenced by the yarn crimp. Hence, the simulated geometry models produced a more accurate prediction of the strain onset of transverse yarn damage accumulation and the stress-strain curve shows the main trend of stiffness reduction. The final failure depends mainly on correct prediction of longitudinal stresses, although the level of damage present at the point of failure affects maximum stress which can be sustained.
Comparison of the two 52% o-FVF models clearly demonstrates the effect of textile geometry on the resulting mechanical Table 3 Comparison of modulus (E) and ultimate tensile strength (S t ) predicted by the models and from experiment.  properties. Non-conservative results were produced by the idealised model because of the geometrical assumptions outlined in Section 2, most notably the lack of yarn waviness. Assessment of the two simulated geometry models, considering the normalised data for the low compaction model, shows a reduction in mechanical properties with compaction since the highly compacted model features greater waviness. In these cases, the relative reduction in tensile strength was over double the relative reduction in stiffness, illustrating that strength is more sensitive to geometrical defects than stiffness. This effect will be even more severe for compressive strength since the presence of waviness causes kink-band formation in yarns [5,6]. Even when considering the raw, un-normalised data of the low o-FVF simulated geometry model, the results suggest that the mechanical properties do not necessarily improve with compaction since the effect of greater o-FVF is counteracted by the reinforcement acting with lower efficiency. This indicates that there is likely to be an optimum moulded thickness for a given reinforcement type and desired properties with regards to a bias towards stiffness or strength. The use of deformation models in conjunction with failure models outlined in this paper could be used help determine this thickness for a given weave style, since the behaviour of 3D woven composites is very dependent on the specific reinforcement and results from one textile cannot be considered to carry over to a different textile. One would expect the fully compacted simulated geometry model to produce the closest match with experiment, however predicted properties were slightly underestimated. One potential source of error was the applied boundary conditions. A unit cell with periodic boundary conditions behaves as if it is part of an infinite material, although due to the physically large unit cell size specimens for testing were only around one unit cell wide for weft loading and two for warp loading. Running unit cell models with no periodicity constraints in the lateral direction however showed negligible changes in stiffness, thus eliminating this as a possible cause. The results of these models are sensitive to material inputs and whilst every care was taken in determining these properties, it is possible that inaccuracies exist between the real material and properties used due to the reliance on manufacturer's and literature data. Fibre transverse modulus for example has a notable effect on overall stiffness and published values vary widely [35,36]. However this was not attributed as the main source of error.
As further investigation, and in order to assess the mesh sensitivity of the models, a convergence study was conducted, generating coarse and refined models by dividing and multiplying the baseline element dimensions by a factor of 1.5 giving element sizes of 0.18 and 0.08 mm respectively. This was done for weft loading in the fully compacted realistic model since there was more damage prior to final failure and this would indicate any differences between different meshes most clearly. The results are summarised in Table 4. The predicted modulus was similar for the coarse and baseline models, but increasing refinement showed a small increase in modulus. This trend of increasing stiffness with voxel refinement occurs for yarns with misalignment to the mesh direction since the voxel discretisation gives a stepped surface which is structurally less efficient. However, when considering straight yarns aligned to the voxel edges, voxel models produce an accurate stiffness prediction compared to a traditional FE approach, even with low mesh refinement. High refinement of misaligned yarns will cause the prediction to converge to the correct numerical solution, though the level of refinement used in this paper is coarser than what is required to achieve full convergence. The implication of this is that the voxel method accurately modelled the idealised  geometry but underestimated stiffness for the two more accurate simulated geometry models. The damage parameter c 2 /c 1 was kept constant in each model such that the strength predictions in the convergence study illustrated the effect of mesh refinement on the damage model. Stiffness degradation due to damage could therefore be expected to retard with refinement because the energy dissipated reduces. However, the effect of refinement on strength did not follow this trend, where in fact there was no clear pattern with both coarse and refined models having a slightly higher strength than the baseline. Extraction of the force through a single weft yarn in the centre of the model at 0.9% strain showed a similar trend to the ultimate strength predictions. The variation in strength was largely due to the variation in voxel discretisation for the different mesh sizes causing variations in the local stress state which fluctuate with mesh refinement rather than converging on a single value. It can thus be concluded that for the range of meshes investigated here, the damage formulation itself is not highly sensitive to the mesh size, once an initial set of parameters is chosen, with voxel discretisation effects dominating. This is a limitation of the voxel method when applied to complex geometrical models, since the key advantage of the technique is the ability to mesh complex geometry. Several techniques could be implemented to improve the accuracy of voxel modelling at low refinement. These include a better representation of the yarn-matrix interface through local refinement [15] or mesh smoothing [37], or by assigning mixed properties to voxels which occupy an interphase region consisting of both yarn and matrix [38,39].

Conclusions
Three different geometrical models of a complex orthogonal 3D woven composite have been assessed using the voxel method and a continuum damage model for determination of elastic moduli and final failure. The assumptions used in an idealised model have been shown to produce significantly non-conservative results for stiffness and especially strength predictions. Models considering realistic geometry from simulation at two levels of compaction produced reduced, and more accurate, predictions. The fully compacted simulated geometry model produced reduced properties compared to the normalised lower compaction simulated geometry model due to the increase in waviness levels. The significant variation of results between the different models clearly illustrates the importance of textile geometry and that idealisation assumptions are not adequate for textiles with complex internal architecture. The fully compacted model should have produced closest match with experimental data, though results were slightly conservative. This has been attributed to the voxel discretisation of the geometry, as the stepped yarn surfaces in the simulated geometry models caused a reduction in stiffness and some fluctuation in strength. Considering such complex geometry introduces greater requirements on mesh refinement in order to achieve convergence than it was feasible to achieve for the considered unit cell and alternative strategies other than mesh refinement need to be considered to mitigate its effects.