Multi-scale modelling of strongly heterogeneous 3D composite structures using spatial Voronoi tessellation

3D composite materials are characterized by complex internal yarn architectures, leading to complex deformation and failure development mechanisms. Net-shaped preforms, which are originally periodic in nature, lose their periodicity when the fabric is draped, deformed on a tool, and consolidated to create geometrically complex composite components. As a result, the internal yarn architecture, which dominates the mechanical behaviour, becomes dependent on the structural geometry. Hence, predicting the mechanical behaviour of 3D composites requires an accurate representation of the yarn architecture within structural scale models. When applied to 3D composites, conventional finite element modelling techniques are limited to either homogenised properties at the structural scale, or the unit cell scale for a more detailed material property definition. Consequently, these models fail to capture the complex phenomena occurring across multiple length scales and their effects on a 3D composite's mechanical response. Here a multi-scale modelling approach based on a 3D spatial Voronoi tessellation is proposed. The model creates an intermediate length scale suitable for homogenization to deal with the non-periodic nature of the final material. Information is passed between the different length scales to allow for the effect of the structural geometry to be taken into account on the smaller scales. The stiffness and surface strain predictions from the proposed model have been found to be in good agreement with experimental results. The proposed modelling framework has been used to gain important insight into the behaviour of this category of materials. It has been observed that the strain and stress distributions are

strongly dependent on the internal yarn architecture and consequently on the final component geometry. Even for simple coupon tests, the internal architecture and geometric effects dominate the mechanical response. Consequently, the behaviour of 3D woven composites should be considered to be a structure specific response rather than generic homogenised material properties.

-Introduction
The increased use of high performance composite materials in a variety of applications requires the adoption of new technologies to improve the performance and reduce the cost of these materials.
The drive for better mechanical performance as well as manufacturing efficiency has led to the introduction of composite materials with through-thickness reinforcement produced in a near netshape pre-form. An example of such a material, is 3D woven composites which are formed from multiple layers of fibre reinforcements woven together with through-thickness binder yarns [1,2].
The result is a woven fabric preform that can be consolidated, infused with resin and cured into the desired final shape. These materials exhibit enhanced impact performance and energy absorption characteristics [3][4][5]. However, the presence of binder yarns during weaving introduces localized deformations such as yarn crimp and waviness [6,7]. Figure [1a, 1b] shows an example of these internal feature on the meso-scale (scale of individual yarns). Additional deformation occurs during the second phase of manufacturing which is the compaction in a mould tool. Unlike the periodic deformations, which are introduced during weaving, the compaction deformation is a result of forming the fibres into the final structural geometric shape. These deformations are thus dependent on the tool geometry and are typically non-periodic in nature. The yarn paths change as the preform is compressed to conform to the tool surface. Additionally, depending on the part final geometry, the preform will experience different levels of compaction at various locations [8]. Figure [1c,1d] shows an example of such deformations on the sub-component or feature scale. Research has shown that the meso-scopic deformations have a significant impact on 3D woven materials mechanical performance [9]; since, the damage initiation process in these materials is dominated by the stress concentrations associated with this localized deformation. Figure [ The yarn architecture deformations occurring during forming and compaction, which will be termed "compaction deformation" hereafter, when combined with the internal features present from material weaving, results in a non-periodic highly heterogeneous material with its properties dependent on the structural geometry. With the help of modelling or experiments the internal yarn architecture of a 3D composite structure can be accurately described. However, this knowledge does not mean that the mechanical behaviour of the final structure can be predicted from finite element models. For conventional solid mechanics problems, it is reasonable to make the separation between the material properties and the structure. The material behaviour is assumed to follow a consistent behaviour throughout certain domains and hence can be represented by an equivalent material model. This kind of separation is not possible for 3D composites since the material behaviour is dependent on the final structural geometry and the manufacturing process.
Additionally, for this category of materials, the strain variation at component scale is comparable to the size of the characteristic model building block. For this case, the material mechanical response cannot be measured consistently using coupon testing at a scale larger than that of the composites constituents, which are the yarn and matrix. On the other hand, building detailed models from the ground up including all the details of the yarn architecture for a full scale structure is not possible with the computational power available today or in the near future. Here, the need arises for a multi-scale modelling approach that can link the basic constituents through multiple length scales all the way to the structural scale. Multi-scale approaches have seen wide use in predicting the behaviour of heterogeneous materials.
Moreover, several of those multi-scale techniques have been adopted to composite materials [10].
These approaches try to address two main challenges arising from the nature of composite materials. The first challenge is including the effect of the yarn architecture in a macro-scale model of the structural behaviour without having to model the heterogeneities explicitly. This is usually achieved by taking a periodic Representative Volume Element (RVE) extracted from the heterogeneous structure. The boundary conditions on RVE's boundaries can be set under an assumption on how this volume interacts with the rest of the material. Periodicity or symmetry of the strain field at the yarn scale are common assumptions. A set of boundary value problems with these conditions provides homogenised/effective properties [11,12]. Single or multiple RVEs can be used to represent the various possible heterogeneities' patterns [13]. This category of approaches has been applied successfully to various types of heterogeneous materials and across multiple length scales. Woven materials introduce an additional level of complexity due to the presence of an intermediate level for the heterogeneity which is the yarn/tow level, here termed the meso-scale.
Dedicated modelling approaches derived from RVE homogenization have been developed to predict the equivalent mechanical properties of a woven unit cell [7,[14][15][16][17]. The underlying assumptions for the RVE approach require that, the heterogeneities or at least the heterogeneities' distribution be periodic [18]. If the structure is irregular or strain field is non periodical, the classical homogenisation framework is not applicable. As shown earlier, for complex geometry 3D woven composite material properties are dependent on both the local tow/yarn architecture and the structural geometry and hence are non-periodic. The results obtained using the RVE approach will thus not be descriptive of the true structural response on the macro-scale.
The other main challenge for multi-scale modelling is the non-linear behaviour of heterogonous materials undergoing progressive damage. Damage initiation at the meso-and micro-scales will progress based on the heterogeneities' patterns and the stress and/or strain state on this scale. The damage will degrade the material on the meso scale, which in turn will impact the global structural behaviour. The interaction between the structure at the macro scale and the damage at the mesoand micro-scales is a nonlinear multi-scale problem. A wide array of multi-scale modelling techniques have been developed for this purpose [19][20][21][22][23][24]. These techniques utilize some form of domain decomposition where the problem domain is divided into two or more domains. The domains can be completely overlapping, partially overlapping or non-overlapping, based on the multi-scale approach used. The smaller scale models can include progressive damage models and are linked to the macro-scale to ensure that the forces and displacements are compatible between the various scales. The modelling approaches in this category are generally efficient and can handle the progressive damage problems well. However, for most these approach the problem domain needs to be divided into areas of interest, which are modelled in detail, while the rest of the domain is homogenized. Decomposing the 3D woven structure before the analysis start with no prior knowledge of the stress state is not a sound strategy. The effect of the local features, which are different throughout the structure, on the mechanical performance cannot be assessed before solving the macro-scale structure. On the other hand, having detailed models for all the possible sub domains will lead to unpractical model sizes where the multi-scale modelling approach losses its efficiency in terms of computational loads.
The topic of damage initiation and progression in composite materials has been widely studied in literature. Physically based phenomenological damage initiation models have been proposed and widely applied for both fibre and matrix dominated properties [25][26][27]. Several continuum damage models based on smeared crack approaches have been proposed for modelling of progressive failure [28][29][30][31][32]. Also, the cohesive zone interface element approach has been used to model delamination and cracking in composite materials [33][34][35][36]. XFEM based methods have also been proposed as an efficient approach to modelling crack growth [37] and have been applied to composite materials [38]. The choice of which model or group of models to use with composite materials is largely related to a balance between accuracy and computational efficiency and is still an open research topic. However, all these models share in common that a detailed knowledge of the stresses and strains is required throughout the domain where the damage model is applied. For 3D woven materials a dedicated multi-scale approach that takes into consideration the complexities associated yarn architecture and its interaction at both the meso and macro scales is needed. The key features for a successful modelling approach, which leads to the understanding of the complex multi-scale nature of 3D woven composites are: • The creation of a new basis of homogenization to replace the absence of the unit cell periodicity.
• Representing the effect of the meso-scopic features on the macro-scale model.
• The ability to locate regions of the problem domain subject to high stress / strain with no prior assumptions as to the whereabouts of these regions within the structure.
• The ability to include detailed meso-scale models with high fidelity at critical locations if required.

-Modelling Approach Overview
The primary goal of a solid mechanics problem is the prediction of mechanical behaviour at the macro-scale where the loads are applied and design criteria can be implemented. Hence, it is essential, when only the properties of the micro scale constituents are available, to link all the length scales from the fibre to the structural levels in the models. In the proposed approach, three length scales will be investigated micro, macro and meso. The micro scale is the scale of the fibre where the material properties are known and where the structure geometry has no significant impact on the mechanical behaviour. The primary driver of the mechanical behaviour at this scale is how tightly the fibres are packed together and how much of the material volume is filled by matrix. The fibre packing is expressed by the Intra Yarn Volume Fraction (IYVF). Micro-mechanical models based on different IYVF can be built using analytical methods or Finite Elements (FE) [39,40]. Based on the IYVF at a given location, a set of equivalent 3D orthotropic material properties can be calculated.
The intermediate scale in this approach is the meso-scale. On the meso-scale, the woven composite structure is described in terms of yarn and matrix materials. The material regions enclosed inside a yarn surface are no longer represented as fibre and matrix but as the homogenized orthotropic equivalent properties, as calculated on the micro-scale for a given IYVF. On this scale, the weave pattern is represented by the yarn paths and yarn cross-sections throughout the fabric. In addition, to the weave pattern and the IYVF, meso-scopic features such as crimp and waviness need to be represented to account for the stress concentration resulting from such features. The third scale is the macro-scale, which is the structure or feature scale. At this scale, the structure geometry is modelled, loads and boundary conditions are applied. Additionally, at the macro-scale, materials are no longer presented as the orthotropic yarn and homogenous matrix but as an equivalent orthotropic material. These macro equivalent material models are calculated from the meso-models for a given sub-domain within the structure. However, since the meso-scopic features are not only a result of the weave pattern but also a result of the structure geometry, the macro and meso scales cannot be separated in the same manner as the separation between the meso and micro-scales.
The multi-scale modelling approach proposed here efficiently handles the modelling of 3D woven structures by solving the problem in two steps. An initial phase where a full scale homogenized mechanical model is generated, taking into account the knowledge of the detailed yarn architecture and the heterogeneities, hereafter called the macro-scale model. The novel macro-scale mechanical model used in this paper will be presented in section 5. The solution of the macro model is used to determine the critical regions within the structure. Next, these critical regions are replaced with high fidelity models where the yarn architecture and heterogeneities are described in sufficient detail, here after called meso-scale models. The meso-scale models are connected to the macro-scale model using a set of Lagrangian Multipliers. This meso/macro-scale interaction allows information to be passed between the two scales thus allowing for the stress state on the meso-scale to affect the response on the macro-scale. Figure

-Internal Yarn Architecture
In the proposed modelling framework, the starting point is the meso-scale models, which connects directly from the geometric modelling phase. The meso-scale models are built based on detailed knowledge of the internal yarn architecture at the given location. These meso models in turn inform the macro models and control the overall structural response. Hence, it is of importance to accurately predict the internal geometry after weaving and compaction.
Several approaches have been adopted to determine the yarn geometries of woven materials. These approaches can be classified, based on the formulation used to describe the dry yarn behaviour, into geometric, kinematic and mechanical models. Geometric approaches assume that the yarn behaviour follows that of geometric entities such as splines [41,42]. A set of geometric rules are used to create the yarn surface geometries ensuring that the yarn paths follow the restriction enforced by the weave pattern and avoiding yarn interpenetration. Such approaches can offer acceptable results for single layer 2D woven materials. However, the coupling of out of plane and in plane deformation in 3D woven fabric create significant challenges for these models that results in less accurate yarn geometries in terms of their path, cross-sections and fibre volume fractions.
The geometries calculated using geometric approaches can be enhanced by eliminating yarn interpenetrations in an additional step using contact algorithms [43]. Another approach to eliminate the mechanical modelling mesh quality issues arising from these models is by using mesh superposition approaches. In these approaches the yarn and matrix material are meshed independently [44,45]. Additionally, experimental data from CT-Scans can be used to augment the geometric approach and considerably enhance the accuracy of such models using direct measurements [46,47]. The key limitation to using CT-Scan approaches is the need to first manufacture physical samples of the material. Another class of models are kinematic models, where contact algorithms are used to simulate the weaving and compaction of 3D materials. In contrast to the experimental/geometric approach, kinematic models can model the internal yarn architecture with minimal experimental inputs. One notable modelling approach is the Digital Element where each yarn is represented by a bundle of beam elements in a contact model [48][49][50].
A digital element model under periodic boundary conditions can be used to accurately simulate the weaving deformations. However, the digital element is normally associated with a high computational effort which limits its applicability to the unit cell scale. Another type of kinematic models is the single surface approach which have been specifically developed to handle large scale compaction models [8,51]. The third category of compaction models is the mechanical models [52][53][54][55]. Several constitutive models for representing the mechanical behaviour of dry fibre yarns have been developed and implemented. These models have the advantage of capturing the compaction pressures and forces as well as the final yarn architecture, which is valuable information from a manufacturing engineering point of view. However, the application of these models are still limited to the unit cell scale due to their computational cost.
For calculating the mechanical performance of 3D woven structures, we are only concerned with determining the yarn/weave architecture as an input to the mechanical models. Consequently, kinematic modelling is the optimal approach for predicting the internal yarn geometry. For the large scale models used in this work, a reduced yarn geometric representation is adopted where each yarn is represented by a single surface. Each yarn is divided into a number of connected segments. Each segment is bound by two yarn sections which are defined as polygons. Segments are straight between the two sections and the yarn surface is linearly interpolated between the two sections. An example of how the fabric yarn architecture is defined in this approach is given in Figure [

-Micro Mechanical Models
Mechanical properties for this level are combined from both micro-scale material models and the yarn architectures.
A key input to the micro-scale models is the IYVF calculated from the meso-level geometry. For each yarn segment as defined in Figure [

-Meso Mechanical Models
Once the equivalent 3D orthotropic material properties have been found for every yarn, the next step is to build a finite element mesh to represent the domain of the problem. Several

-Full scale models using spatial Voronoi tessellation
The meso-model described in the previous section is not practical from a computational expense point of view to apply to the feature or structural scale. In order to accurately describe the yarn surface geometries and resin pockets, a high fidelity mesh with a large number of elements is required. To avoid this cost we propose a macro scale modelling approach using a Voronoi tessellation as a basis for homogenization. Voronoi tessellation is a numerical algorithm to divide a spatial domain into completely interlocking cells which tessellate to form the original domain. The conventional tessellation creates its cells such that for all the points inside a cell with a given centre, the distance to the cell centre is smaller than the distance to any other centre in the given domain. Voronoi tessellation has been used in multi-scale modelling of heterogeneous materials with cells built around the heterogeneities as geometric centres [64][65][66]. Each Voronoi cell is then homogenised taking into account the heterogeneity in the centre. Several challenges face applying the conventional Voronoi tessellation to the 3D composites. First, the yarns follow complex paths which may or may not be related to a specific layer.
Additionally, the yarn cross-section changes for each yarn segment. Hence, it is required to expand the tessellation to account for the 3D nature of the problem and for the variability of the yarn cross-section.
For this work, the conventional 2D point based tessellation has been replaced by a 3D volume based tessellation. In this tessellation, the yarn segment surfaces have been used as tessellation generators. As a result, the cells are the volumes closest to a given segment of a yarn surface. Since each yarn segment is straight, the typical Voronoi cell consists of a single straight yarn segment surrounded by an arbitrary region of matrix material. The precise shape of the cell outer surface is dependent on the yarn segment cross-section shape and also the yarn segments geometry in its vicinity and is generally irregular.
Consequently, the tessellation decomposes the 3D composite structure into a set of interlocking cells each containing a 1D yarn segment. Each cell material properties can be homogenized using volume integration over a very fine grid to calculate the cell homogenized material properties with no need for orientation averaging. The averaged material properties are considered to be acting in same direction as the 1D yarn segment at its centre. These 3D orthotropic material properties are then assigned to a macro voxel mesh of the complete structure. The result of the complete process is that the structure has been

Figure 9. Voronoi tessellation generators a) Yarn segment defined by points b) Yarn segments represented as polyhedron.
For the homogenization to be successful it is essential to have a Voronoi tessellation which satisfies two conditions: • The tessellations cover the complete domain with no regions left without being assigned to a Voronoi cell.
• The Voronoi cells are mutually exclusive which ensures that there is no overlap of the cells to avoid inaccurate homogenization.  Using the proposed tessellation, a problem domain containing any 3D yarn architecture can be completely divided into cells which serves as basis for homogenization. The pseudo-code for the implementation of this homogenization algorithm is shown in Figure [10]. The core of this homogenization algorithm is organizing the faces defining each yarn surface into a search tree. The search tree is used to calculate the homogenized properties and to assign these properties to the macroscale Voxel mesh. The search trees in this work was implemented using the AABB Tree from the CGAL library. Figure 11. shows the tessellation process applied to a simple reinforcement geometry. The geometry shown here is of two out of phase spiral yarns embedded in a cuboidal domain. The tessellation decomposes the domain into a set of interlocking cells which are assigned to specific yarn sections. Multiple cross-section are taken along the domain length to show how the domain is divided between the two yarns. The same approach used to decompose the two yarns problem in Figure 11 can be used to decompose any 3D woven composite, regardless of its complexity or the number of yarns.

Figure 11. 3D spatial Voronoi tessellation applied to a simple geometry, a) Reinforcements geometry, b) Domain decomposition using Voronoi tessellation, c) Reinforcements' domains after decomposition, d) Cross-sections across the tessellation domain.
For purpose of model verification and validation, a set of meso and macro models have been created and compared against experimental results. Three sets of models have been solved where each set models a specific physical experiment. The models include two tensile samples, one loaded in the warp direction and the other in the weft direction, and a V-notch rail shear test [68]. The high fidelity models, loading directions and the Voronoi models are shown in Figure [12]. For the shear test, the shear modulus was calculated following the ASTM standard D 7078 which uses the cross-section area between the two notches to calculate modulus. The same standard was applied to calculate the shear modulus from the meso and macro scale models. The elastic moduli predicted by both the meso and the macro scales have been compared against the value measured by experiments. The results are shown in Table [2]. Model sizes data are shown in Table [

6-Meso-scale boundary value problems
Once the macro-scale problem has been solved and the displacement field calculated, the highly loaded regions within that model can be determined with reasonable accuracy. At this stage the overall stiffness Where ‫ܭ‬ is the stiffness matrix of the meso-scale model and ߣ is the Lagrangian multiplier forces.
Lagrangian multipliers by definition introduce singularities to the diagonal of the system matrix, which in turn causes most linear equations solvers to fail. As a result, the complete system has to be solved in an iterative process where the Lagrangian multiplier forces are fed back to the macro-scale system. The Lagrangian forces can be projected on the macro-scale mesh using the Moore-Penrose pseudoinverse of the coefficient matrix ܰ ା after scaling for macro/meso mesh size ratio. In this approach, the macroscale degrees of freedom can be divided into internal and boundary degrees of freedom following the FETI approach [19,69]. Then the macro system can be rewritten as: Here, ‫ܭ‬ ெ is the stiffness matrix of the internal degrees of freedom. ‫ܭ‬ ெ is the stiffness matrix components associated with the degrees of freedom on the boundary between the meso and macro models. ‫ܭ‬ ெ and ‫ܭ‬ ெ is the macro stiffness boundary/internal coupling terms. ‫ݑ‬ ெ is the displacement associated with the internal degrees of freedom on the macro-scale. ‫ݑ‬ ெ is the displacement associated with the boundary degrees of freedom. ݂ ெ is applied forces associated with the macro-scale internal degrees of freedom and ݂ ெ is applied forces associated with meso/macro boundary. Figure [19]. The comparison shows that there is no detectable displacement discontinuity between the macro and meso regions in the hybrid model. A key drawback of many global-local analysis approaches is the presence of discontinuities between the global and local scales. These discontinuities can lead to stress artefacts in the meso-models. For the proposed modelling approach, since the macro displacement field is heterogeneous and follows the same pattern as the meso scale, no such stress artefacts can be seen in the meso model. Figure [20] shows displacement field on the interface between the meso and macro scale ‫ݑ(‬ ெ ) which has been calculated by solving equations [2] and [3] simultaneously. Figure [21] shows the stress distribution on the middle slice as predicted by both the multi-scale and the full meso scale model. Figure Figure [24b]. Figure [25] shows the multiscale mesh and the displacement as calculated by the multi-scale model. The full macro-scale model is 2 million degrees of freedom and the hybrid meso/macro is 3.2 million degrees of freedom. Figure [ 26] displays the fibre direction stresses in the meso-model. Figure [27

-Conclusions
In this paper, a novel integrated multi-scale modelling framework for 3D woven structures has been proposed and verified against experimental results. The proposed model has been used to gain insight into the mechanical behaviour of 3D woven composites. The strain and stress distributions found using this modelling approach show a strong relation between the internal yarn architecture and the mechanical response of 3D woven composites. Internal yarn deformations such as crimp and waviness create localized stress concentrations that can act as damage initiators in these structures. Additionally, it has been shown that the internal yarn architecture changes throughout the manufacturing phases and is dependent on the final structural geometry. Consequently, the yarn architecture and the meso-scale mechanical response become no longer only material properties but are also influenced by the specific structure of interest. This observation leads to additional complications when modelling or experimentally testing 3D composites due to the need to consider the yarn architecture at the mesoscale, even for structural scale response.
Composites modelling have traditionally relied on achieving scale separation between the macro and meso scale. Here, it has been shown that for 3D composites the loss of periodicity and the mechanical response dependency on yarn architecture mean that scale separation is no longer feasible. A coupled approach which links the meso-scale to the global structural response is required. This coupled approach needs to include the influence of the internal yarn architecture on both the macro and meso scales.
Additionally, coupon testing is a widely used approach of composites characterization. In this work, we have modelled three different specimens commonly used in such tests. The strain profiles calculated by the models show considerable variation throughout the specimens' planform and through thickness.
Even for simple test such as a tensile test, the strain distribution shows large variation across the sample width, which is an indicator of strong interaction between the yarn architecture and the sample size.
These yarn architecture dependent variations mean that conventional coupon testing of these materials will measure specimen specific properties rather than generic material properties. For 3D woven composite structures with complex geometry, the changes in internal architecture are associated with a changes to the yarn architecture on the meso-scale. As a results, homogenised material properties measured from flat samples or calculated from models of flat samples will not be valid as inputs to structural scale models. This effect has been observed clearly in the double curvature component modelled in this paper. In this model, the highest transverse stresses were observed in regions where the internal yarn architecture showed the highest deformations, as a result of the component geometry.
This observation shows that more elaborate specimen design, which is informed by the internal architecture and the final structure geometry, is needed to carry out proper characterisation. The modelling techniques presented here offer an opportunity to carry out trials of such characterisations in a virtual sense, due to their accurate physical representation of geometry at the relevant scales.