(2018). An iterative multiscale modelling approach for nonlinear analysis of 3D composites.

The advent of new more complex classes of strongly heterogeneous materials, such as 3D woven com- posites, introduces new challenges for well-established ﬁnite element multiscale modelling approaches. These materials’ internal architecture dominates the local stress concentrations, damage initiation and damage progression. Additionally, the material loses periodicity during manufacture and conventional homogenization approaches become inapplicable. In this paper, a multiscale modelling approach based on domain decomposition and homogenization is proposed to model the mechanical behaviour of these materials. The proposed model formulates a set of displacement and force compatibility conditions between the various subdomains. The compatibility conditions are formulated by limiting the set of kinematically admissible solutions on the smaller scales in order to satisfy the larger scale basis functions at the inter- faces. The different subdomains are solved alongside the compatibility conditions in an iterative process. The proposed multiscale framework reduces the stress artefacts on the subdomains’ boundaries. This fea- ture allows the 3D woven material internal architecture represented at the smaller scales to control the structural response at the global scale. Additionally, this framework allows for selective application of nonlinear material models to the subdomains of interest through its ability to redistribute stresses across the subdomains boundaries.


Introduction
3D woven composites are manufactured by integrally weaving multiple 2D layers of fibre reinforcement together with fibre bundles (yarns) in the through-thickness direction. The complete woven preform is then compacted to the desired thickness and shape. Next, the preforms are infused using a liquid resin matrix material, which is then cured under pressure. Compared to 2D composites, 3D composites offer enhanced impact resistance, delamination resistance and energy absorption characteristics ( Guénon et al., 1989 ). However, the complex nature of the fibre architecture in these composite materials poses significant challenges to conventional modelling approaches. The internal fibre architecture of a complex 3D composite structure can be non-periodic and exhibits extensive features and deformations ( El Said et al., 2014 ). Here features are defined as unavoidable yarn path deformation in order to conform to the geometry and/or the weave architecture.
On the other hand, defects are defined as additional deformation resulting from the manufacturing processes and/or handling such as yarn crimp and waviness. Fig. 1 shows a cross-section of a kinematic simulation of 3D composite forming, demonstrating the loss of periodicity at the structural scale. Research has shown that damage initiation occurs at areas of geometric deformation in the fibre architecture such as crimp and waviness ( Tan et al., 20 0 0 ;Mahadik and Hallett, 2011 ). Moreover, it has been observed that stress relief mechanisms, such as matrix cracks, act to redistribute the stresses and impact damage progression ( Cox et al., 1994( Cox et al., , 1996. Because of these mechanisms, failure strains, stresses and damage tolerance in 3D composites can only be predicted by modelling the interaction of multiple failure events and the material internal yarn architecture.
The process of progressive failure in a 3D woven composite material typically takes place on three different length scales. The first scale is the micro-scale, which is the scale of the individual fibres and the surrounding matrix regions. The intermediate scale is usually termed the meso-scale, which is at the scale of the yarns. The meso-scale is where the weave style is represented and yarn architectural features are described. Finally, the   macro-scale is where the global structure geometry is defined, loads are applied and the design criteria such as structural stiffness can be observed. Fig. 2 defines these length scales for a 3D woven composite material. A large body of published research is dedicated to damage modelling in woven composite materials. A prerequisite to the application of these models is a detailed knowledge of the stress/strain state across the domain where damage models are applied. For 3D woven composites, this stress/strain knowledge cannot be generated using high fidelity finite element models alone, since the loading and damage processes occur across multiple length scales. On the other hand, modelling all the relevant scales in detail would result in prohibitively large models.
In most conventional composites applications, the three length scales can be conveniently separated using periodic homogenization. Equivalent microscale material properties can be calculated based on the knowledge of the fibre properties, matrix properties and the fibre density in a given Representative Volume Element (RVE). This can be done using either analytical models ( Chamis, 1983( Chamis, , 1984 or computational micro-mechanics models ( Arteiro et al., 2015 ;Melro et al., 2013 ;Bai et al., 2015 ). In these RVE models, clusters of individual fibres embedded in a matrix are modelled under periodic boundary conditions. The clustering of these fibres is chosen to represent the fibre density at a given location. The microscale equivalent material properties can then be used to build a mesoscale RVE which is used to calculate another set of material properties taking into account the fibre/yarn orientations and stacking sequences ( Lee et al., 2005 ;Koumpias et al., 2014 ;Obert et al., 2014 ;Schneider et al. ;Green et al., 2014 b;Biragoni and Hallett ;Tsukrov et al., 2011 ). Finally, a macroscale structural model can be built using the mesoscale equivalent properties. This macroscale model takes into consideration the structure geometry, the applied loads and boundary conditions. However, in 3D woven materials, the loss of periodicity on the macroscale means that scale separation cannot be attained between the meso and macro scales. Consequently, more involved multi-scale modelling techniques are needed to handle 3D woven structures.
The topic of multiscale modelling of composites has been widely covered in literature ( Kanouté et al., 2009 ). These models can range from simple ones concerned with calculating envelopes for the elastic properties of composites ( Hill, 1963 ;Hashin and Shtrikman, 1963 ;Watt, 1979 ) to complex numerical models ( Ladevèze and Lubineau, 2001 ;Ladevèze et al., 2006 ). Nonlinear multiscale models can be divided into two main categories. The first category is homogenization based approaches, where meso-scale models are used to calculate equivalent material properties, which can be later used as the constitutive law in macro-scale models. Homogenisation based models are intuitive and straightforward to implement since the constitutive laws are normally built on the mesoscale, independent from the macroscale simulation. Composite materials pose a challenge for homogenisation based models. Progressive damage in composites is a non-linear phenomenon that is dependent on the stress state.
The key challenge in this case is establishing a link between the macro-scale and meso-scale models. This link should allow the macro-stress state to drive the damage on the meso-scale. The second category of models is based on domain decomposition techniques. These techniques divide the model into interconnected regions. Each region or domain can be solved under a different modelling framework. The main computational benefits come from the ability to run the expensive nonlinear material models on the regions of interest rather than the whole structure.
The question of which multiscale approach to choose is dependent on whether the scale separation hypothesis is reasonable. In those cases where scale separation is possible, homogenization approaches are normally the first choice. Several homogenizationbased approaches have been proposed for multiscale modelling of composite materials mechanical problems ( Kouznetsova et al., 2002 ;Kami ński and Kleiber, 2000 ;Chung and Tamma, 1999 ). In the case of linear material behaviour, an RVE can be solved separately from the macro simulation. However, for the case of nonlinear material behaviour, which is the case for damage in composites, the mesoscale response becomes dependent on the macroscale stress state. In such cases, each macroscale integration point can be assigned a mesoscale RVE under periodic boundary conditions as is the case with FE 2 approaches ( Feyel and Chaboche, 20 0 0 ;Feyel, 2003 ;Smit et al., 1998 ). The stresses and forces on the macro-scale are calculated as the volume average of the stress and forces on the RVE. The boundary conditions on the RVE are found from the macro-scale solution and applied to the meso-scale. The feedback loop between the meso and macro scale is maintained with each FE solver iteration until the simulation is complete. These approaches are quite powerful since they are nonlinear in the true sense. However, the need for an RVE per integration point results in a computationally expensive simulation. Additionally, the numerical errors in these methods increase as the RVE size approaches the macro mesh size, which is the case with 3D woven composites. The smallest periodic unit in a 3D material is normally few centimetres in size, which approaches the size of the geometric features of most advanced engineering components.
Another main category of multi-scale models that are applicable to composites are the global-local analysis approaches. These approaches can be used to model specific subdomains of a given structures at the meso-scale. Finite Element Tearing and Interconnecting (FETI) is an important example of such methods Roux, 1991 , 1992 ;Farhat et al., 1998 , 20 0 0 , 20 01 ). In these methods, the problem domain is decomposed into multiple disconnected subdomains. The subdomains are then re-connected using a set of Lagrangian multipliers introduced to enforce compatibility conditions at the interfaces. The Lagrangian equations and the subdomains' FE problems can be solved iteratively until the sub domains are in equilibrium under the Lagrangian forces. Additionally, each sub domain can be meshed separately with different non-conformal meshes that leads to further runtime reductions. However, the uncoupling of the subdomains in these approaches is based on the decomposition of the subdomains'stiffness matrices. Consequently, while these approaches can handle linear and geometrically nonlinear problems, material non-linearity can be challenging to implement efficiently.
A different approach to non-linear global/local analysis is to use an iterative approach to achieve displacement and force compatibility at the subdomains' interfaces ( Gendre et al., 2009( Gendre et al., , 2011Whitcomb, 1991 ). Several global/local analysis approaches have been applied to woven composites ( Woo and Whitcomb, 1996 ).
A key strength of this approach is the ability to start from a linear solution of the complete structure, which enhances conver-  gence speed. Additionally, these approaches can handle material nonlinearity resulting from damage. However, some of these approaches are limited to subdomains with compatible meshes, which is a serious limitation when modelling materials with complex internal architectures. Other approaches decompose the domain into a set of interlocking sub domains and interfaces. For each interface between two sub domains a displacement compatibility condition and stress compatibility conditions are applied ( Ladevèze et al., , 2010. The two compatibility conditions are achieved simultaneously during the solution process. The meso-scale subdomain displacements and forces are related to the macro quantities on the interface using projector functions. Similar domain decomposition approaches have been applied to crack propagation in complex structures ( Guidault et al.,20 07 ,20 08 ) and composite materials ( Daghia and Ladevèze, 2012 ;Kerfriden et al., 2009 ). An alternative approach is to employ model     reduction techniques and global/local analysis to couple a reduced macroscale model with a mesoscale model undergoing damage ( Kerfriden et al., 2012 ). Thus, focusing the computational cost on the mesoscale region.
A limitation for some of the domain decomposition approaches is the need to impose either a displacement or a force compatibility condition, this limitation can lead to stress artefacts on the sub domain boundaries. For 3D woven composites and other strongly heterogeneous materials, damage initiation starts from architectural features distributed throughout the material. Boundary stress artefacts will distort the damage initiations and progression, and will result in inaccurate predictions. For 3D woven and other similarly complex heterogeneous materials a dedicated multi-scale modelling approach is therefore needed, which  takes into consideration the complexities associated with 3D meso-structure. The proposed model needs to minimize stress artefacts resulting from subdomain connections. Additionally, the proposed framework should be capable of handling the exceptionally large finite element models resulting from the discretization of 3D composite geometries on the meso-scale. Here we start from the non-intrusive iterative global/local analysis approach proposed by Gendre et al. (2009 ) and the iterative approach developed by Daghia and Ladevèze (2012 ) . These approaches will be expanded to account for non-conformal meshes using displacement compatibility conditions based on the macro-scale basis functions. Additionally, in order to reduce meso/macro interface stress artefacts a homogenization technique based on Voronoi tessellation will be used to create the macro-scale linear models . The next sections in this paper are dedicated to presenting this approach and its application to 3D woven composites.

Compatibility conditions
The proposed approach starts by modelling the entire structure with a homogenised formulation, which will be termed the macro-model hereafter. The homogenisation approach used in this work follows the Voronoi tessellation approach proposed by El Said et al. (2016 ). In this approach, a volume based Voronoi tessellation is used to smear the matrix regions in a 3D woven composite while maintaining the local material orientation. The resulting homogenised model is fully anisotropic and accurately represents the structure's stiffness response. This homogenized model is solved under the global loads and boundary conditions. Once the macroscale problem has been solved and the displacement and strain fields calculated, the highly loaded regions within that model can  be determined with reasonable accuracy. At this stage, the overall stiffness response of the component has been calculated and the next stage of analysis is to understand the failure envelope of the structure under investigation. Thus, it becomes necessary to remodel the highly loaded subdomains using a high fidelity representation of the heterogeneities and internal material architecture. These subdomains will replace the equivalent regions in the macro-model.
At first glance, this approach appears to be decomposing the domain in the spatial sense. This spatial division is the philosophy behind the global/local analysis approaches reviewed earlier. However, the efficiency of the proposed approach can be enhanced by augmenting the modelling approach with key elements from homogenization techniques. This is achieved by using a homogenised model to solve the elastic loading phase. Additionally, the homogenised solution is also used as the initial condition for the high-fidelity model in the global/local analysis. Besides the spatial decomposition, the system response will be divided into homogenized and high fidelity solutions. For the subdomains modelled as macro-scale, the proposed approach will be primarily concerned with calculating stiffness and displacement response. The external loads and boundary conditions will be applied on the macroscale model. In addition, at this scale the structural design criteria can be assessed. For the mesoscale subdomains, the material internal architecture will be described in details and the stress/strain state observed. Additionally, at this fine scale the material failure can be predicted. Fig. 3 shows an overview of the proposed multi-scale approach.
The main component of the proposed modelling approach is the connection between the meso and macro scale subdomains. Consider a heterogeneous body under a given set of loads and boundary conditions acting on its surface. The domain can be divided into two sub domains 1 and 2 separated by a common boundary ∂ 12 ( Fig. 4 ). If the sub domain 1 is meshed using a coarse mesh and the subdomain 2 is meshed in a fine mesh, 2 can be solved as mesoscale while 1 can be solved as macroscale (homogenised). By applying this decomposition, the stress and strain distribution in the 2 domain can be predicted at the scale  of individual yarns and damage models applied. At the same time, the stiffness response of 1 is represented in the meso-model. The incompatibility of displacements / forces on the boundaries can lead to artificial stress concentrations at the interface. For 3D woven composites and other strongly heterogeneous materials, the stress concentrations resulting from the internal architecture are distributed across the problem domain. Moreover, significant stress redistribution is associated with the progressive damage development. Consequently, the interfacial stresses will dominate the damage initiation and progression.
Purely descending approaches will not capture the effect of the meso-model during the macro-model computation since the two models are solved sequentially and only once. In a descending approach, one quantity, either displacement or force, is imposed from the macro-model onto the meso-model. The dual quantity (for example, the force in the case of a displacement formulation) will not be continuous. In order to make these quantities continuous, a dual condition is needed which can be written as a residual. This dual condition can then be solved in an iterative process, thus satisfying both conditions simultaneously ensuring force and displacement continuity. To avoid these errors, a combined displacement / force formulation for the connection between the two scales is proposed in this work. This proposed approach employs an iterative scheme to relieve the artificial stresses resulting from the direct coupling of the two domains.
In the case of non-conforming meshes, the number of degrees of freedom on the meso/macro interface do not match. Consequently, a choice has to be made about which quantities should be continuous locally everywhere and which quantities should be continuous only on average. The work by Guidault et al. (2007 ) has investigated the issue of errors arising from mesh discontinuities for displacement and force based formulations. In particular, imposing the macro displacement completely might result in over constraining the meso degrees of freedom on the interface. From this point of view, the choice of completely connecting the meso nodes to the macro displacement, which is the option retained in the present work, is not optimal. However, this approach offers a pragmatic solution to the problem of domain decomposition in 3D woven composites where the internal architecture is continuous across the domain interfaces.
In the proposed approach, two compatibility conditions are assembled and enforced on the sub domains interfaces. First, a displacement compatibility condition is formulated, which establishes a relation between the displacements at the macro and meso scale interfaces. The proposed displacement condition is assembled based on the macro subdomain ( 1 ) finite elements' shape functions. This approach is applicable as long as the following conditions are satisfied: where d 1 and d 2 are the average element size in the domains 1 and 2 , respectively. Consequently, 1 will be considered as the macro domain as long as it is represented by a coarser mesh, regardless of the actual domain size. Next, the set of kinematically admissible solutions of subdomain 2 will be limited to the set of solutions that fit the solution of 1 on the boundary ∂ 12 . The two subdomains are connected by enforcing the macroscale shape functions at the interface between the two subdomains.. Examples of non-conforming meshes on the interface are shown in Fig. 5 . For the displacement of node n 1 to be constrained to the macro element edge N 3 − N 4 , n 1 displacement has to satisfy the macroscale element shape function. The displacement of the meso-scale node n 1 can be interpolated as follows: Here A 1 ( η) to A 4 ( η) are the shape functions associated with each macro node in the quadrilateral macro element, u n 1 is the nodal displacement of the relevant meso-node and U Ni are the macro-nodes' displacements. The shape functions are computed at the point n 1 ( η), which corresponds to the meso-node coordinates described in the macro-element intrinsic coordinates system. For a general 3D macro element with M nodes, the meso-node displacement can be written as: Alternatively, for all nodes on the interface ∂ 12 : where A is a l × M coefficient matrix for an interface with l meso-nodes and M macro-nodes. The compatibility condition provided by Eq. (2) establishes a relation between the meso and macro displacement. Since, this compatibility condition links meso and macro nodes it is independent from the macro/meso elements' relations. Fig. 5 shows three different arrangements between macro and meso models. Each of these arrangements can be handled by the proposed approach. Fig. 6 shows a schematic of feasible and unfeasible meso-scale boundary conditions based on the kinematics of the macro mesh.
A similar compatibility condition is needed to connect the macro and meso interfacial forces acting between the two subdomains. This condition can be developed by enforcing a force balance at the sub domains' interface using projection of the mesoscale forces. In this approach, the meso-scale nodal forces are projected on the macro interface elements using, once again, the macro-scale shape functions. For instance, the contribution of the meso-scale nodal force from node n 1 in Fig. 6 to the macro-scale nodal forces is given by: where A 1 ( η) to A 4 ( η) have been previously defined, f n 1 is the meso-scale force applied to n 1 and F n 1 Ni are the contributions of f n 1 projected on the macro-scale nodal forces F Ni . Summation of the contributions of all meso-scale nodal forces yields the relation between the macro-scale and meso-scale nodal forces over the whole interface ∂ 12 : where F ∂ 12 is the global force vector acting on the macro-scale interface, f ∂ 12 is the global force vector acting on meso-scale interface and A T is the transpose of the coefficient matrix A described for the displacement boundary conditions. The meso-scale nodal force on the boundary is calculated by integrating the stresses on the interfacial elements in similar manner to assembling an internal force vector.
Additionally, it can be shown that the proposed projection satisfies the work equivalence at the macro/meso interface. The work on the meso-domain boundary is given by: Similarly, the macro domain boundary is given by: Once the two conditions are assembled, Eqs. (2) and (3) couple the sub-domains without selecting a force based or displacement based formulation. Hence, the problem has been divided into two sub domains, with associated compatibility conditions.

Computational strategy
The main goal of the proposed computational strategy is to solve the multiple domains under given loads and boundary conditions, whilst satisfying the compatibility conditions given by Eqs.
(2) and (3) . The complete system is non-linear and cannot be solved directly using linear equation solvers. Here we propose an iterative scheme to solve this non-linear system. An initial homogenised model of the complete problem is built. This homogenised model is then solved to locate the highly stressed areas in the structure. Once these highly stressed regions are located, these regions are remodelled using the mesoscale formulation and a fine mesh. Next, the compatibility conditions are assembled for the macro and meso-scale models. Then, the macro-scale shape functions are used to assemble the coefficients matrix A . The meso-scale subdomain is solved as boundary value problem under the displacements calculated from the full macro-scale solution using Eq. (2) . Next, a new set of interfacial forces from the mesomodel can be extracted and projected on the macro boundary using Eq. (3) . A residual force vector for iteration n can be calculated by comparing the change in the macro/meso interface forces: Fig. 7 shows a flow chart of this computational strategy. The iterative process between the meso and macro scales can be repeated until the residual force vector converges. To achieve convergence, Eq. (3) is not applied directly in the numerical implementation. Instead, an interpolation scheme is used to calculate the interfacial forces, taking into account the forces from current iteration ( n ) and the previous iteration ( n − 1) using a factor α, and is given by: An iterative sparse conjugate gradient solver from the PETSC package ( Balay et al., 20016 , 2016 ) is used to solve the finite element equations. The convergence speed for this category of solvers is dependent on the initial condition assumed at the start of the iterations. Here, we can project the homogenised solution from the first iteration on the mesoscale domain and use it as an initial condition for the iterative solver. The projection can be done using the same kinematic admissibility conditions, which were used to assemble the displacement compatibility conditions described in Eq. (2) . However, here the projection will be applied for the complete meso-scale domain: Here the term δu 1 represents the mesoscale solution complement with respect to the homogenized solution. Using this approach, the multiscale problem will be about finding the difference between the full homogenised solution and the multiscale model rather than solving the complete model from scratch, which leads to faster convergence. A key strength of this approach is avoiding the need to calculate the inverse or pseudo-inverse of the subdomain stiffness matrix. Consequently, it is not required to assemble a global stiffness matrix in the case of explicit time integration to solve the multiscale problem.
A quarter model of square homogenous plate with an open hole at its centre is used as a benchmark problem to illustrate the iterative strategy. The plate, under tension, is modelled in displacement control. The domain decomposition and the solution progress are shown in Figs. 8 and 9 . The region around the hole is modelled using a fine mesh while the rest of the plate is modelled  using a coarse mesh. The two meshes are incompatible with a size ratio of 1:6. Figs. 10 and 11 show the evolution of the average displacement and stresses calculated for every cross-section along the plate length. It can be seen that as the solution progress both the stresses and displacement become continuous across the plate. Thus, eliminating the stress jump that existed at the start of the solution. Fig. 12 shows a comparison of the stresses between the multiscale solution and a full mesoscale solution, where the whole plate is modelled using a fine mesh. The maximum error in the mesoscale region is around 3%. This is reasonable value given the benefits of eliminating stress discontinuity, which will be demonstrated when applying the models to 3D woven composites.

Time dependent problems
The domain decomposition approach described in the previous section can be expanded to handle time dependent problems. This can be achieved by updating the displacement compatibility conditions to link the velocities on the interface of the meso and macro models. Two different time integration schemes can be used for this category of problems: implicit and explicit. An implicit time integration scheme requires the assembly of a Jacobian matrix for each of the subdomains (meso and macro) and achieving convergence for each iteration. Under the current strategy, the Jacobian assembly and the residual vector are calculated independently for each domain. Convergence for the time integration scheme is also checked independently for each domain. Once both domains achieve convergence, the compatibility conditions can be updated and the multi-scale convergence is checked. Fig. 13 shows the flow chart for the time integration schemes for the proposed multi-scale solver.
Progressive damage modelling is an essential aspect of 3D woven composite modelling. The damage initiation stresses in these materials can be quite low. Due to the heterogeneous nature of composites and its ability to redistribute stresses, these materials will continue to carry the loads well after damage has initiated. One of the challenges for progressive damage finite element models in composites is the handling of fully damaged elements using continuum damage models. As damage progresses through the model, it is common for some elements to be completely damaged before a final failure event. These completely failed elements can lead to instability issues in an implicit time integration scheme. It is therefore desirable, when modelling progressive failure in 3D woven composites, to be able to switch to explicit time integration when instabilities are detected. The time step in an explicit integration is considerably smaller than implicit schemes. Consequently, applying the multiscale convergence criteria for every time step will result in unreasonable runtimes. For explicit time integration, the compatibility conditions are updated based on velocities and force values from the previous time step. Consequently, the multiscale convergence check can be omitted in this case. Fig. 13 also shows the integration scheme for explicit time integration. A dedicated parallel finite element code has been developed to implement the multiscale approach described here including the automatic switching from implicit to explicit time integration as needed.

Application to 3D woven composites
The modelling framework described so far is generic and can be applied to any non-linear finite element problem. In this section, the model will be applied to 3D woven composites to model their mechanical performance, based on realistic internal architectures. A prerequisite to modelling heterogeneous materials is having a complete understanding of the heterogeneities' distributions. This can be achieved in 3D composites by employing kinematic modelling, which simulates the weaving and compaction processes. In this paper, the kinematic modelling approach developed by El Said et al. (2014 ) was used to model a 3D orthogonal woven composite described by Green et al. (2014 a). Fig. 14 shows the internal yarn architecture of the 3D woven fabric used in this paper, as predicted by a multi-scale kinematic model. Once the internal material architecture is known, the material orientations, fibre volume fractions and matrix material regions can be determined and applied to a finite element model. Here, a Voxel meshing approach is used to handle the geometric complexities associated with 3D woven architectures. The construction of these high fidelity finite element models using Voxelization has been discussed at length in Green et al. (2014 b) andEl Said et al. (2016 ). For 3D woven composites with complex yarn architectures, the number of features and their geometric complexity becomes prohibitive for most meshing techniques. Voxelization divides the full problem domain into equally sized cubes called Voxels. Each individual Voxel is compared to the geometric entities in the domain to define its properties. Using this approach any complex geometry can be discretized across the domain. In the case of two phase materials, such as 3D woven composites, the yarns can be mapped to the Voxels, with the remaining Voxels assigned the matrix material properties. An example of the Voxelization of a 3D woven composite is shown in Fig. 14 (b) and (c). The presence of spurious stress concentrations on the matrix/ yarn boundaries is sometimes associated with Voxelization. In this paper, a direct mapping approach is used to assign material properties to each integration point. In this approach, Voxels which are modeled as fully integrated Hexahedral elements, can have matrix and yarn properties at the same time. This approach reduces spurious stress concentrations by creating transition elements on the surface of each yarn . It is worth mentioning that other techniques have been proposed in literature to handle spurious Voxelization stress concentrations, including non-local stress averaging and damage tracking ( Fang et al., 2016 ).
In this work, yarn elements will be assigned 3D orthotropic material properties oriented in the direction of the fibres at each location. By choosing to assign equivalent material properties to the Voxels, an implicit assumption of scale separation, between the microscale (fibre scale) and mesoscale (yarn scale), is made. In this case, the assigned properties have to be representative of the microscale fibre and matrix properties. Several homogenised models have been proposed in literature and seen wide use for this purpose. In this work, the model by Chamis (1983Chamis ( , 1984 is employed to find microscale homogenised properties taking into account the microscale fibre volume fraction for each Voxel. In later sections of this paper, a homogenized damage model is introduced for non-linear analysis. This is also done under an assumption of micro/meso scale separation. There is however no requirement to have full micro / meso scale separation in the proposed framework. Nonlinear multiscale modelling approaches, such as the ones described in the introduction of this paper, can also be used to link the micro and meso scales.
The voxel models described for the mesoscale are capable of predicting the stress concentrations resulting from architectural features. This is achieved by describing the yarn paths, crimp, waviness and matrix regions with a high fidelity. Associated with this high fidelity description is a high computational cost for these models, which means they cannot be used on the macroscale even for linear elastic analysis. For the 3D woven examples in this paper, a spatial surface based Voronoi tessellation proposed by El  is used to construct the homogenized macroscale models. Voronoi tessellation homogenization has seen wide use in modelling multi-phase materials by decomposing the domain into cells centred on the inclusions. Equivalent properties are calculated for each cell accounting for the inclusion and surrounding matrix. The surface based Voronoi tessellation has the added advantage of considering the inclusion morphology in the tessellation. Thus, the homogenised model maintains the material directionality and the fibre density. This is of prime importance in the case of global/local analysis of 3D woven composites. Using this homogenization approach, the material directions in the macro subdomains will match those across the interface in the meso subdomains, allowing the macro model to respond more flexibly to the changes in the meso models. This leads to fewer stress artefacts and better stress redistribution across the boundary that will be demonstrated in Section 4 of this paper. Fig. 15 shows the volume based Voronoi tessellation applied to a 3D composites sample.
The framework described earlier in this section, which consists of a macroscale Voronoi tessellation model, a Voxel based mesoscale model and the iterative global / local analysis, has been applied to several 3D woven composites examples. Fig. 16 shows the domain decomposition of a 3D composite V-notch shear specimen ( Adams et al., 2007 ). This specimen is under in-plane shear loading and the critical section is between the notch tips. The critical section surrounding the notch has been modelled as meso-scale while the rest of the problem is modelled as macro-scale. A comparison between the resultant strain predicted by the multi-scale model and a full high fidelity solution is also shown in Fig. 17 . It can be seen that for the multi-scale model no discontinuities are present at the boundaries between the meso and macro scales. Moreover, the stiffness response of the two models matches well, with the shear modulus predicted by the multi-scale model being 4.7 GPa and by the full mesoscale model 4.8 GPa.

Multi-scale modelling of progressive damage
At this stage, a detailed prediction of the stresses and strains across regions of interest in a 3D woven composite can be generated. This can be used to further predict the progressive damage behaviour of the material. 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 ( Puck and Schürmann, 1998 ;Cuntze and Freund, 2004 ;Caddell et al., 1974 ). Several continuum damage models based on smeared crack approaches have been proposed for modelling of progressive failure ( Pinho et al., 2006b, Donadon et al., 2008Vogler et al., 2013 ;Camanho et al., 2013 ;Lapczyk and Hurtado, 2007 ). In addition, cohesive zone models have been used to model delamination and cracking ( Jiang et al., 2007 ;De Moura and Gonçalves, 2004 ;Mukhopadhyay et al., 2015 a, b ). XFEM based methods have also been proposed as an efficient approach to modelling crack growth ( Moës and Belytschko, 2002 ) and have been applied to composite materials ( Ye et al., 2012 ). 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, it is essential to understand the different length scales that comprise a given multiscale problem and determine the scale at which the damage or the nonlinearity can be introduced. In this paper, scale separation has been assumed between the meso and micro scales. Consequently, the microscale behaviour can be represented in the mesoscale model through equivalent properties, which includes parameters such as strength and energy release rates.
In the case of 3D woven composites, two damage models are needed at the meso-scale. A damage model representing the yarn material phase, which at the microscale is composed of fibres embedded in a matrix material. The other damage model represents the matrix regions, which are composed of homogeneous material. The selection of suitable scale to represent nonlinearity is specific to each multiscale problem. It is not mandatory to apply the material nonlinearity to the meso-scale. For other types of materials, it could be necessary to implement the damage models on finer scales or it might be sufficient to apply them on the macro-scale. In the next section, the implementation of damage models to the mesoscale, based on equivalent microscale properties, is presented. as this is suitable for the 3D woven composites under consideration here. These damage models were selected to meet the requirements for the verification case presented in this paper, which is an in-plane tension loading case. However, for different loading cases and a more robust general solution, more complex failure criteria will be needed. For an instance, to model compression loading cases, the compressive fibre failure criterion should include the effect of shear stress, such as in Pinho et al. (2006 a).

Energy regularized smeared crack models
A smeared crack damage model is chosen for use in both the yarn and matrix materials, since these models are computationally efficient and suitable for Voxel meshes. A main drawback of smeared crack models is that the mesh can affect the direction of damage progression. The cubic Voxel meshes used here tend to reduce these mesh bias effects. Additionally, theses damage models assume an element crack plane. The cubic Voxels make it easier to define this plane. It should be noted that these models can be applied to other types of 3D elements if required. Smeared crack models are straightforward to implement for the mesoscale, where the yarns are represented as continuum entities. In this paper, the smeared crack model proposed by Pinho et al. (2006 b) and expanded to composites containing ply waviness by Mukhopadhyay et al. (2015a , b ) is used. A brief description of this model is given next for purpose of completeness. The interested reader is referred to the original references for more details. This model tracks the initiation and progression of the various damage modes in a unidirectional composite, which is representative of the yarn materials in 3D woven composites. The model assumes that damage initiates in the form of a crack on the microscale, which later progresses to split a mesoscale element completely.
The key damage modes' initiation and progression are defined in relation to the stress components acting on an assumed crack surface. Fig. 18 shows the definition of these stresses in relation to the crack surface and crack angle. σ n is the normal stress perpendicular to crack surface. τ T and τ L are the in-plane shear stress acting on the crack surface. The angle ϕ is the angle between the crack surface and the material axis.
Undamaged matrix response for yarn and resin pockets: in this model the matrix regions for both yarn and resin pockets material models is assumed to have a non-linear out of plane shearing behaviour governed by the following relation ( Mukhopadhyay et al., 2015 b): Here σ 1 i is the shear stress components, γ 1 i is the associated shear strain, A and B are material constants. Eq. (6) describes a nonlinear shear stress/strain relationship during loading. For the unloading regime, the material behaviour is assumed linear.
Fibre failure: a tensile criterion based on fibre rupture for the yarn material is given by: where S c 11 is the mesoscale fibre ultimate tensile strength, which has been calculated from the microscale volume fraction at the given location following Chamis, 1983Chamis, , 1984 ) based on the material properties of the constituents given in Green et al. (2014 b). The finite element implementation of this criterion reduces the fibre direction stresses after failure by a factor of 0.98. The maximum reduction factor has been limited to 0.98 to avoid numerical instabilities.  8) and (9) are calculated on the fracture plane, which is not known beforehand. A golden section search algorithm ( Fang et al., 2016 ) is used to find the fracture angle ϕ which maximizes the failure index.

Matrix crack initiation under compression (in yarn material): matrix crack initiation in compression is
Matrix crack progressive damage model using smeared crack approach: Once damage initiation for matrix cracking is predicted, an energy regularized smeared crack model is used to predict the crack opening. After the crack initiation, a driving stress for the crack propagation is calculated: Similarly, a driving strain ε mat and a crack surface shear strain γ mat can be calculated by resolving the strain components in relation to the crack surface. Using those values, an effective mixed mode critical fracture energy per unit area is calculated at the point of crack initiation: where G m is the mixed mode critical fracture energy per unit area, G N is the fracture energy per unit area in the normal direction acting perpendicular to the fracture plane, G T and G L are the shear fracture energies per unit area acting in the fracture plane. σ 0 and ε 0 are the effective stress and strains at the point of crack initiation acting in the relevant direction. Assuming a bilinear damage propagation law on the meso-scale, an ultimate failure strain can be calculated from the area under the curve in Fig. 19 , which is given by: The characteristic length l c ensures mesh independent energy dissipation. This length is an element property, calculated by projecting the length of the element side on the crack surface. A damage index can be calculated relating the stress and strain on the crack surface is defined by: The stresses acting on the crack surface can be adjusted to account for the progressive damage: To avoid numerical instability, in the finite element implementation, the value of D mc is limited to 0.98.
Resin pockets -matrix material progressive damage model using smeared crack: the damage model discussed so far covers the yarn transverse cracking and fibre failure for yarn materials, which are orthotropic. In 3D woven composites, the mechanical performance is dominated by the yarn behaviour, due to the presence of yarns acting in all three directions. The mechanical performance of the matrix regions is thus of reduced importance as compared to 2D woven composites. In this work, a simplified damage model for the resin pockets is employed. This model uses a pressure dependent Von Mises criteria following Green et al. (2014 b) for damage initiation ( Eq. (15 )).
where σ 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. The influence of hydrostatic pressure in the tensile case evaluated here is however negligible.
Damage progression is controlled by the same smeared crack model used yarn materials. The fracture plane in the matrix domain is oriented using two angles instead of one and the golden section search has to be adapted accordingly. Using this search, the plane maximizing the matrix principal stresses is determined. Once the fracture plane is found, the smeared crack model can then be applied in a manner similar to the yarn material matrix crack.

Damage model results and analysis
The smeared damage model described in the previous section has been introduced to the proposed iterative multiscale approach. A test model of an open hole 3D woven composite tension specimen was run. The specimen is 29.7 mm in width and 5.3 mm in thickness. The 3D woven composite used in this specimen is a carbon fibre orthogonal 5 harness satin weave. The details of this fabric architecture and material properties is given in El Said et al. (2014 and Green et al. (2014 a). The domain surrounding the open hole was modelled as mesoscale while the rest of the specimen was modelled as macroscale. Fig. 20 shows the domain decomposition into meso and macroscale. The damage model was applied to the mesoscale domain only while the macro domain was homogenized using Voronoi tessellation. The model was loaded in tension under displacement control.
The combined multi-scale/progressive damage capability is an effective tool to understand the mechanical behaviour of complex 3D composite materials. Fig. 21 shows the mesoscale section binder yarn ( z -direction) transverse matrix cracks at increasing load levels. It can be seen that the transverse cracks appear at the regions with highest level of internal architectural defects (crimp). The initiation takes place around 30% Ultimate Tensile Strength (UTS). Damage initiation at this stage appears to be driven by the internal fibre architecture only, independent of the location of the hole. At higher loading levels, it can be seen that damage progresses faster in regions where the binder yarns interact with the geometric features (the open hole). In domain decomposition approaches, where stress artefacts appear on the subdomain's interfaces, the damage may initiate from the boundaries instead of from the material's meso-scale features. Consequently, the interaction with geometric features and the ultimate failure behaviour will not be predicted correctly. Figs. 22 and 23 show the evolution of average stress and displacement across the model. The results show no significant stress or displacement discontinuities across the meso/macro boundaries. Fig. 24 shows post-mortem images of two physical specimens of an open hole tension test. The images show a highly diffused damage pattern in the samples. Additionally, the images show a high variability in damage pattern between the samples. This variability is also reflected in the UTS values measured experimentally (76-94 kN) ( Townsend et al., 2009 ). Fig. 25 shows a comparison of the experimentally measured resultant forces of several 3D woven specimens and the simulation which further confirms the variability. This observation can be attributed to the interaction between the material internal architecture and the geometric features mentioned earlier. The damage pattern predicted by the model for matrix cracking, in both yarns and resin pockets, shows a similarly diffuse pattern. Additionally, the simulated damage patterns replicate the features observed in the experiment. Surface cracks and yarn splitting originating from the hole was observed in both the model and experiments over matching paths.
The multi-scale framework captures another important aspect of 3D woven composites' mechanical behaviour, which is their ability to redistribute stresses away from damaged zones. Fig. 26 shows the fibre direction stress distribution in the open hole tension model across the meso and macro sub domains. Initially, stress concentrations appear near the hole, which leads to first fibre failure at this location. However, this does not lead to catastrophic failure straight away. The material manages to divert the stresses from the broken fibres, loading other undamaged fibres, which allows the material to continue carrying load. During this phase, the iterative approach allows the surface yarns in the macro subdomain to react to the increase of stresses in the surface yarns in the meso subdomain. Eventually, the level of stress in the fibres leads to their complete rupture across the specimen and ultimate failure. The multi-scale models presented here were capable of capturing the ultimate failure load with good accuracy, given the level of variability existing in these materials. Key to this accurate prediction is the successful redistribution of stresses across the subdomain boundaries, which is a result of the iterative approach proposed here.

Discussion and conclusions
In this paper, a multi-scale modelling approach for the prediction of progressive damage in 3D woven materials is presented. It starts by solving a homogenized model of the complete structure. Next, the critically loaded regions in the structure are remodelled using a mesoscale high fidelity formulation, where damage models can be applied. The meso/macro problem is then solved using an iterative multi-scale process. The iterative process employs displacement and force compatibility conditions to connect the various subdomains and relieve the boundary stress artefacts. The compatibility conditions are formulated by limiting the set of kinematically admissible solutions on the smaller scales in order to satisfy the larger scale basis functions at the interfaces.
The proposed modelling approach was combined with a progressive damage model to predict the failure behaviour of 3D woven composites. It has been demonstrated that the proposed modelling framework is a capable of handling the strongly heterogeneous behaviour of these types of materials. The stress concentrations from the mesoscale material architecture were predicted and allowed to initiate damage, free from boundary stress artefacts. Additionally, the interaction between the material architecture and the geometric features in the structure has been captured. Moreover, the proposed model tracked the stress redistribution in the material undergoing progressive damage across the different subdomains due interactions between the problem subdomains throughout the various iterations. While this approach has been demonstrated here in relation to 3D woven composites in a meso/macro model, it is equally applicable to other categories of strongly heterogeneous materials.
The model results and experimental verification demonstrate the non-periodic nature of 3D woven composite structures. The meso-scale stress concentrations interact with the geometric features and lead to a loss periodicity in the stress fields. Traditionally, composites have been modelled using periodic homogenisation. The work presented here indicates that periodic homogenisation of 3D woven composites is not so valuable at the macro-scale as for other composite architectures. To illustrate this, the V-notch specimen model presented earlier has been repeated using material properties obtained from periodic homogenisation and is shown in Fig. 27 . A unit cell model under periodic boundary conditions was used to calculate the equivalent properties of the 3D woven material. Six load cases were solved to calculate the equivalent orthotropic material properties. These properties were then applied to a macro-scale model and solved using ABAQUS standard. It can be seen that the homogenised model fails to capture the stress concentration resulting from the internal architecture and its interaction with the geometric features. While these types of homogenised models can predict the stiffness behaviour with reasonable accuracy, they fail to predict the local variations in stress state and consequently are not adequate for damage and failure modelling.