Meso-scale modelling of 3D woven composite T-joints with weave variations

: A meso-scale modelling framework is proposed to simulate the 3D woven fibre architectures and the mechanical performance of the composite T-joints, subjected to quasi-static tensile pull-off loading. The proposed method starts with building the realistic reinforcement geometries of the 3D woven T-joints at the mesoscale, of which the modelling strategy is applicable for other types of geometries with weave variations at the T-joint junction. Damage modelling incorporates both interface and constituent material damage, in conjunction with a continuum damage mechanics approach to account for the progressive failure behaviour. With a voxel based cohesive zone model, the proposed method is able to model mode I delamination based on the voxel mesh technique, which has advantages in meshing. Predicted results are in good agreement with experimental data beyond initial failure, in terms of load-displacement responses, failure events, damage initiation and propagation. The significant effect of fibre architecture variations on mechanical behaviour is successfully predicted through this modelling method without any further correlation of input parameters in damage model. This predictive method will facilitate the design and optimisation of 3D woven T-joint preforms.


Introduction
For 3D woven composite structures, especially for those with geometric features, the design space of their preforms is large with an enormous amount of variations in the 3D spatial reinforcement architecture.Understanding the influence of the fibre architecture of 3D woven composites on their mechanical properties is fundamental to the design phase.However, at present this is mainly dependent on experimental testing [1][2][3][4], due to the lack of analysis techniques that are able to predict the resulting structural performance for 3D weave architectures, which restricts the application of 3D woven composites.
The hierarchy of textile composites is usually classified according to the length scales: fibres

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
2 in matrix of the impregnated yarns at the micro-scale, impregnated textile reinforcements and bulk matrix at the meso-scale and 3D composite components at the macro-scale.Usually the impregnated yarns are locally considered as unidirectional (UD) composites so that most of the theories for modelling of UD composites are still applicable for textile composites at the meso-scale.Multiscale modelling techniques have therefore been widely used in modelling textile composite structures.Among the steps summarised by Lomov et al. [5], the meshing difficulty persists, due to complex fibre architecture in textile composites.Geometry simplification, such as artificially reducing the size of yarn cross-sections to eliminate extreme thin layers of matrix in-between adjacent yarns [5][6][7], is usually employed but would lead to the usage of a higher intra-yarn fibre volume fraction as well as an unrealistic constituent interface.Alternatively, voxel-based FE method has been proved to be an effective way in stress/strain analysis of textile composites for their significant advantage in meshing [8][9][10], albeit spurious prediction on damage initiation for a multi-layer plain woven composite was found by [11].Mesh dependency was found by Ernst et al. [10] when using conformal mesh to analyse the failure of textile composites with fracture energy approach, due to elements with irregular aspect ratios would be usually generated near the constituent interface if the interface is not formed of flat surfaces as seen in textile composites.Instead, voxel mesh was adopted as the mesh dependency vanished and also good agreement in stiffness and progressive damage analysis was observed between simulations and experimental results for a thick NCF specimen subjected to three-point bending load [10].
Apparently, voxel method has both its cons and pros and sometimes a compromise of using voxel mesh has to be taken for modelling of composites with complex fibre architecture when conformal mesh is not readily available based on the state-of-art meshing technique whilst voxel mesh is capable to achieve most of the required result.
For damage in directions other than the fibre direction, brittle failure does not always occur

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
3 and thus continuum damage mechanics (CDM) are often exploited to model damage propagation.Damage initiation is first evaluated by a failure criteria and then a degradation scheme on elastic constants is applied to the stiffness matrix.For instance, Nobeen et al. [12] proposed a modelling method for the progressive damage of braided composites based on an instantaneous damage model, in which an instantaneous reduction (degradation factor) on the material constants was applied once the damage for the impregnated yarns was predicted by Hashin's criteria.It should be noted that when damage was predicted in the fibre direction, the constants were reduced to a near zero value as it was assumed to be a complete failure.Good agreement with experiments was observed, however, one limitation with this method is that the degradation factor needs to be correlated with the experimental results.Similar schemes to degrade the elastic constants was also previously used in [13,14].Puck and Schurmann [15] developed a phenomenological law to degrade the material constants after damage initiation, in conjunction with their failure criterion.Progressive damage behaviour can be predicted by this method and they offered recommendations on the selection of empirical parameters in the damage model in the absence of experimental data.A similar phenomenological damage model was also proposed by Ruijter [16].Although there are suggestions on selection of empirical parameters in these models, they still need to be fully validated against experimental data due to the phenomenological nature.Energy-based damage models were proposed in several studies [17,18], but this requires the determination of fracture toughness through experiment as an input to the models.
A number of meso-scale FE models based on simple flat unit cells for the mechanical performance of 3D woven composites showing good agreement with experiments were reported [19][20][21], but most of them were not being validated for a different weave pattern to justify the predictive capability.For 3D woven composites with geometric features, the fibre architecture would become more complex than for 3D woven composite flat panels, and no publications on the modelling of such materials at meso-scale have been reported.This work is to develop a meso-scale modelling method that can predict the mechanical behaviour of different 3D woven composite T-joints due to weave variation subjected to quasi-static tensile pull-off loading.The meso-scale T-joint model that reflects the feature of the reinforcement architecture is built based on the geometric modelling strategy introduced in Section 3. Details including boundary conditions and constituent properties along with damage modelling techniques are given in Sections 4-6.The results are compared with experimental data in terms of load-displacement response, failure modes, damage initiation and propagation in Section 7.

Materials and testing
Two types of 3D woven T-joint preform were used in this study manufactured by Sigmatex based on Hexcel IM7 12K carbon fibre.The preforms are based on a 3D orthogonal weave with the only variation at the junction.Specimens were woven flat with pre-positioned bifurcations on a Jacquard machine and then folded into a T shape.clamped at two ends by the fixture [22].A single-lens DANTEC Q400 Digital Image Correlation (DIC) system was used to monitor the full-field strains around the junction regions of all specimens.The obtained testing results [22] will be used to validate the proposed modelling method in Section 7.

Construction of meso-scale 3D woven T-joint models
Generating the reinforcement geometry within a composite plays a pivotal part in the mesoscale FE analysis.The accuracy of the predicted composite performance can be improved based on a meso-scale model with realistic fibre geometry [8,23].The reinforcement geometries of the two types of 3D woven T-joint were modelled using TexGen [24], based on the geometric parameters of yarns extracted from µCT analysis [25].The extracted CT measurements with a resolution of 30 µm/pixel include yarn cross-section dimensions, yarn spacing, cross-section centre points location along yarn paths.The fibre architectures of the two types of 3D woven T-joint have identical 3D orthogonal weave patterns in the flange and web sections, with the only difference being the geometry of the junction region.Thus construction of the meso-scale models for the composite T-joints followed the strategy of dividing the T-shaped structure into three sub-geometries: the junction region representing  Benefiting from the periodicity in woven reinforcements, only a repeat unit (unit cell) was modelled along the width direction (6.7mm in z-axis in Fig. 2), which is about one third of the specimen width (20 mm).Additionally, the length of the web in the geometric models is reduced, as the deformation in the web was negligible compared with the deflection from bending of the flange, due to the high effective modulus in the web along its length/loading direction (y-axis in Fig. 2).A voxel meshing method was used in discretising the geometric models due to the complexity in the fibre architectures of the 3D woven T-joints where conformal meshing is difficult to achieve.Fig. 3 shows the voxel-discretised models generated by TexGen based on the above geometries.They were discretised with an element length of 0.1 mm after a mesh convergence study on a reduced model for elasticity and damage (Fig. 6).Cohesive surfaces between yarns and matrix were then added in the voxel models using Hypermesh.Note that in order to reduce the computational cost, there are no interface elements added (perfectly bonded interface) in the web of the 3D woven type 2 model because no delamination was observed in the web region during the experimental tests.

Boundary conditions
To simulate the quasi-static tensile pull-off tests on the 3D woven T-joint specimens, a displacement-controlled load (smooth amplitude) in the y-axis direction was prescribed on the nodes of the top surface of the models, with the other two translational degree of freedom constrained (Fig. 4).AA'DD' and BB'CC' are the surfaces at the edges of clamps on the flange of the specimens, and therefore the nodes on them were assumed to be fixed.To avoid over-constraint on the nodes of the boundaries of the flange, the interfacial nodes of the yarns (nodes on yarn cross-section boundaries) on surfaces AA'DD' and BB'CC' were not constrained in the above manner, and thus they were free in the boundary conditions but constrained by the cohesive tractions from the corresponding interfacial nodes of the matrix.
Periodic boundary conditions are expected to be prescribed on the front and back surfaces of the unit cell models.But this was shown to lead to a higher computation cost because of the large size of the FE models.In addition, a deviation of less than 2% in elastic response between unit cell models with and without the periodic boundary conditions was found through a study of braided composites under uniaxial loading [6].Similarly, non-periodic boundary conditions were also used by other authors in the failure analysis of textile composites [12], as the requirement of identical coordinates for the node pairs on the opposite faces of the model for periodic boundary conditions is difficult to meet for textile composites.
In this study, the effects of periodic boundary conditions on structural stiffness and strength was studied on a simplified model for the 3D woven T-joint (with weft yarns only) subjected to the same loading condition [26].The results for a unit cell width model with and without periodic BCs in the width direction were compared to the results of full-size model.A deviation of about 5% in stiffness and strength was identified and therefore periodic boundary conditions were not adopted in the 3D woven T-joint models.

Determination of intra-yarn fibre volume fraction
An averaged intra-yarn fibre volume fraction (intra-yarn V f ) by preserving the overall V f of the composite based on the volume of yarns in the specific geometric model is commonly used to characterise the yarn's properties in modelling of textile composites [27].However, the intra-yarn V f may vary in different yarns or at different locations of a yarn.Instead of using an averaged intra-yarn V f for the whole model, variation in the intra-yarn V f for warp, weft and binder yarns were considered, which were respectively calculated by matching the V f of the composite in each yarn direction.An approximation of the fibre volume fraction of the composite in the weft yarn direction ܸ ௪௧ is as follows, assuming a straight yarn path (without crimp): where ‫ܣ‬ ௦௧ is the area of the cross-section of the flange; ‫ܣ‬ ௬ = 12݇ × ߨ݀ ଶ /4 is the total fibre area in the section for a weft yarn of 12K filament count, ݀ is the fibre diameter; ݊ ௪௧ is the total number of weft yarns in the flange which can be obtained from the µCT scan (Fig. 7).Then the intra-yarn V f for the weft yarns can be obtained through dividing the  The method above is also applicable for the warp yarns.A coefficient greater than one is needed on the numerator of Eq. ( 1) to account for the yarn path curvature of the binders.µCT scan is not always necessary if ݊ ௪௧ can instead be acquired from the preform manufacturer.
Based on Eq. ( 1), for the 3D woven type 1 T-joint model, the calculated intra-yarn V f for the weft yarns is 62.7%, which is close to that of the warp (71.9%) and binder yarns (56.8%).To simplify the model, the same intra-yarn V f of the weft yarns was used across the whole model in the FE analysis.Because the predicted failure behaviour is far more sensitive to the properties of the weft yarns than those of the warp or binder yarns, a small deviation is likely to be caused by ignoring the variation in the intra-yarn fibre volume fraction.At the same time, the intra-yarn V f was considered unchanged for the two types of 3D woven T-joint because they were made of the same materials at the same global fibre volume fraction.

Homogenized yarn properties
The properties of the yarn elements are approximated by the properties of a homogenized UD composite with a V f equivalent to the intra-yarn V f .Similar to dealing with laminates, the homogenized yarn is assumed to have transversely isotropic properties which can be obtained either by the micro-scale unit cell modelling or an analytical solution based on micromechanics.The Chamis model [29] was used to calculated the homogenized yarn properties.The calculated properties along with the properties for Hexcel IM7 carbon fibre and Gurit Prime 20LV epoxy resin are listed in Table 1.Due to the absence of experimental data, a set of empirical formulae for calculating the strengths of the homogenized yarn was used [5]: where F is the material strength and superscripts t and c denote tension and compression.
Table 2 lists the calculated strengths of the homogenised yarn for using in the FE analysis.

Damage modelling
To model the failure of composites, the damage modes considered in the FE models should cover all the potential failure modes, although sometimes for simplicity, only the failure events observed in mechanical tests are included.In this study, the damage modelling incorporated both yarn/matrix interface damage and damage in constituent materials, i.e. bulk matrix and homogenized yarn materials.

Interface damage modelling
Delamination was found to be a typical failure mode in the testing of composite T-joints and cohesive zone model (CZM) been proven to be effective to model the delamination [30,31].
However, it is open to doubt if using surface-based CZM with voxel method to model delamination is feasible, as the interface elements will be generated on the step-like constituent boundaries.Zhang el al. [9] used voxel mesh to perform damage simulation of a braided composite with surface-based CZM accounting for tow-tow delamination and good correlation with experimental results was obtained.However, some authors stated that it is not possible to use voxel mesh in CZM, as generation of interface elements on step-like interface would be problematic [32], or the interface damage initiation and fracture energy cannot be computed on a step-like interface [33].Before proceeding with computationally intensive studies on the 3D woven T-joint models, a reduced T-joint model, with each half comprising of four layers of bent uniaxial non-crimp fabric without any fixation material, was used to study the effects of meshing technique on the mechanical performance of the structure with CZM.Both conformal mesh (hexahedral C3D8R) and voxel mesh (C3D8R) were generated for this geometry as shown in Fig. 6.The boundary conditions, load case and material properties on the reduce model are assumed to be the same as the 3D woven T-joint models.(0.08 mm).Therefore 0.1 mm was adopted as baseline element lengths for the following voxel-based CZM analysis.In addition, the stress contours of the two models at a same displacement were compared and good agreement was also observed.
The surface-based cohesive behaviour in Abaqus was used in this study, formulated by the bilinear constitutive law.It is noted that other more accurate constitutive laws [34,35] for crack propagation were proposed recently, but due to commercial availability, the quadratic stress criterion and mixed mode power law were selected for damage initiation and evolution: where ߬ is the initial failure stress, ‫ܩ‬ is critical energy release rate for each of the single delamination modes, with subscript stands for normal, first and second shear directions; and p is the power in the criterion.'˂˃' is the Macaulay operator.
The interface stiffness (slope for the elastic stage of the bilinear law) could affect the global compliance prior to damage initiation.As it is not a material constant but a numerical value introduced by the CZM method, the selection of interface stiffness value was initially empirical or based on correlation with experimental results.An analytical approach for estimating the minimum/converged interface stiffness applicable for textile composites under model I delamination was proposed, by extending a previous study for laminates only [36] (given in Supplementary data).In addition, the effect of voxel mesh interface on structural stiffness was analysed in terms of elastic loading and damage initiation according to the interface formulation (Supplementary data).from the elastic loading stage was quantified but it is difficult to evaluate the contribution from premature interface damage analytically.
Numerical investigation was then performed based on the reduced structure introduced before.Interface properties used were similar with previous studies for carbon/epoxy composites [30,31] but with varied interface stiffness values (Fig. 7).Compared with the results from the equivalent conformal model, the voxel model was able to capture the failure load but overestimated the structural stiffness, which agrees with the analysis in supplementary data.Through the parametric study on interface stiffness, it was found that if the interface stiffness was reduced by one order of magnitude (0.1) from the converged magnitude (10 5 MPa/mm) to compensate the overestimation in structural stiffness, the voxel model showed a similar load-displacement response to the conformal model.Therefore these properties (Table 3) will be used in the 3D woven T-joint modelling and validated against experiment data in the following section.It should be noted that interfacial strength and toughness also affect the load-displacement response of the T-joint, but they are limited to the failure stage and would not the change to initial stiffness of the response.In addition, this paper used the interface properties characterised from specimens with a flat interface MPa/mm 30 MPa 60 MPa 0.22 mJ/mm 2 1.2 mJ/mm 2 1.2 mJ/mm 2 2

Constituent material damage modelling
For damage in the homogenized yarn material, Hashin's failure criteria [37] developed for UD composites were used here to capture the damage initiation for each failure mode: ; ; ; where I 1 to I 4 are the four stress invariants; Bulk matrix damage was evaluated by the pressure dependent modified von Mises criterion which can take into account the difference between tensile and compressive strength for the isotropic material [38]: where ‫ܨ‬ , ‫ܨ‬ ௧ are the compressive and tensile strengths of the bulk matrix respectively.
After damage initiation, the behaviour of the damaged constituent materials was modelled by a CDM approach, which degrades the moduli of the damaged constituent materials through a phenomenological law first proposed by Ruijter [16]: where P(d i ) is a stiffness penalty factor function for degrading the corresponding modulus in terms of the failure modes d i defined in Eq. ( 5) and ( 6).c 1 and c 2 are constants and Ruijter

M A N U S C R I P T A C C E P T E D
ACCEPTED MANUSCRIPT [16] found that c 1 =8 and c 2 =13 gave good agreement with the experimental stress-strain response for plain weave carbon/epoxy composites under tensile load.A minimum value of 0.001 for the stiffness penalty factor P(d i ) was maintained when the material is considered to be fully damaged to avoid numerical instability.
Thus the elastic properties of the damaged yarn material are: where E, G with superscript d denote the moduli of the damaged yarn material.Note that catastrophic damage is assumed after the initiation of compressive failure (d 2 ≥1) in the fibre direction, as in Hashin's failure criteria (Eq.( 5)) it is only determined by the stress invariant I 1 .The properties relating to the transverse direction are also gradually degraded with the accumulation of damage in the fibre direction.Poisson's ratios ν 12 , ν 13 and ν 23 are assumed to be unchanged in order to maintain a symmetric stiffness matrix after damage initiation.
Similarly, the Young's modulus of the damaged bulk matrix material is: where E with superscript d denotes the Young's modulus of the damaged matrix material.
The stiffness penalty factor for the Young's modulus of damaged bulk matrix is calculated by using d m and the same constants c 1 and c 2 in Eq. (7).Mesh dependency should be avoided when modelling composites damage with CDM [39].This was assured by comparing the results with a finer mesh (0.08 mm) based on the reduced T-joint geometry introduced in Section 6.1.

Results and discussion
FE models were solved by Abaqus Explicit 6.13 with a user-defined material subroutine (VUMAT) for matrix and yarn materials with behaviour defined in Section 6.2.The predicted behaviour from the proposed FE method for modelling the two types of 3D woven T-joint is properties were used in the FE models for simulation of the tensile pull-off tests.Fig. 8 shows the comparison of predicted and experimental load-displacement responses for the two Tjoints respectively.After introducing the voxel-based cohesive surface, the non-linearity in the stiffness was captured accurately, whilst a previous study based on perfect-bonding condition was found to over-predict the stiffness [40].It should be noted that it is difficult to determine an accurate loading displacement for damage initiation unless from the load-displacement responses, as damage initiated inside the specimens is hard to observe in the testing and also that the first failure of an element (material or interface) in the FE models could not represent the macroscopic damage in the specimens.In the tests, the damage onset of the type 2 specimen occurred later than the type 1 specimen leading to a higher initial failure load, and this feature resulting from the weave variation was successfully captured by the FE method.From the analysis of the test results, it   Besides the accurate prediction capability for the behaviour in the elastic and part of the progressive damage stages, the power of the proposed modelling method is that it does not require any model input parameter correlation with experimental results when the fibre architecture is changed at the junction.High efficiency from this voxel-based method should also be highlighted when modelling composites with complex fibre reinforcements, where conformal meshing is difficult to achieve.

Conclusions
A meso-scale modelling method to predict weave architectures of 3D woven composite Tjoints and the resulting mechanical behaviour under quasi-static tensile pull-off loading was proposed, results were shown to agree well with experimental data beyond initial failure.The proposed method starts with a strategy building the meso-scale geometries of the T-joints, of which the modelling strategy is applicable for other types of geometries with weave variations at the junction.Damage modelling incorporated both yarn/matrix interface damage and damage in bulk matrix and homogenized yarn materials, in conjunction with a continuum damage mechanics approach to account for the progressive failure behaviour.Predicted results are in good agreement with experimental data beyond initial failure, in terms of loaddisplacement responses, failure events, damage initiation and propagation.However, premature catastrophic failure for the two models was found in elements at flange boundary and this is likely to be caused by the simplification of the boundary conditions.

Fig. 1 (
(a), (b)) from x-ray micro computed tomography (µCT) shows the fibre architectures with the direction of weft yarns marked, illustrating 3D woven type 2, where half of the weft yarns are crossing over the other half at the junction, in comparison with 3D woven type 1.The two types of 3D woven composite T-joint specimen reinforced by the above preforms were moulded through a vacuum assisted resin transfer moulding process infused with Gurit Prime 20LV epoxy resin, giving a fibre volume fraction of 45%.The T-joint specimens were cut and tested under quasi-static tensile pull-off loading (Fig. 1(c) and (d)).The clamps were a custom-designed fixture in stainless steel with M6 bolts.There are about 20 mm length of T-joint flange clamped into each side.This fixture was subsequently bolted onto a steel I-Beam attached to an Instron 5581 test machine with a static 50 kN load cell.A displacement load at a rate of 1mm/min was applied on the web of the specimen (three tests for each type), with the flange M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 5

Fig. 1 .
Fig. 1.Images from µCT scan of the two types of 3D preforms showing the weave variation at the junction: (a) Type 1; (b) Type 2; (c) testing layout for T-joints; (d) specimen geometry and fixture

6 the
weave variation, and the flange and web sections which are two unit cells of 3D orthogonal weaves, as illustrated in Fig. 2 (left) for the geometry of the 3D woven type 1 Tjoint.This modelling strategy facilitates the construction of other types of 3D woven T-joints by only varying the geometry of the junction region (Fig. 2 (right)).For Type 2 joint, the weft yarns go into out-of-page direction due to weave variation.The CT image (Figure 1 b) was a 2D slicing, while the geometric model (Figure 2 left) was a 3D projection view.

Fig. 4 .
Fig. 4. BCs for the 3D woven T-joint models subjected to a quasi-static tensile pull-off load M A N U S C R I P T A C C E P T E D above fibre volume fraction of the composite in the weft yarn direction ܸ ௪௧ by the volume fraction of weft yarns in the voxel model.

Fig. 5 .
Fig. 5. µCT scan of the section of the flange showing ݊ ௪௧ =30, there are 6 yarn stacks and 5 weft yarns per stack

Fig. 6 .
Fig. 6.Reduced model geometry and its conformal (element length: approx.0.08) and solid mesh (element length: 0.1), units: mmA mesh convergence study, simplified as an elastic analysis under perfectly bonded The analysis was limited to model I delamination as this is the critical case in T-joint loading.It is found that voxel discretisation would change the stress state at the interface and consequently the voxel model would show a stiffer behaviour than the conformal model.The additional stiffness for the voxel model resulting M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 13

d 1 and d 2
are the damage parameters for fibre dominated failure modes; d 3 and d 4 are the damage parameters for transverse matrix dominated failure modes.
E P T E D compared to the experimental results, in terms of initial stiffness, damage imitation, and failure mode and damage propagation.Because both T-joint specimens are made of same materials along with the same fibre volume fraction, identical material and interface

Fig. 8 .
Fig. 8. Predicted load-displacement responses for the type 1 (left) and type 2 (right) specimens in comparison with test results, FE stopped due to convergence problems (enlarged views of initial section shown on the bottom row)

Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .
Fig. 9. Comparison of progressive interface damage in 3D woven type 1 (left) and 3D woven type 2 (right) FE models; the failed interfaces (full failure) are shown in red As shown in Fig. 10, the delamination onset locations in the type 1 model coincide with those The proposed method is able to model mode I delamination based on the voxel mesh technique, which is advantageous in meshing.More significantly, it does not require any parameter correlation in damage model with experimental results when the fibre architecture is changed at the M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 21 junction, new weave patterns of the 3D woven T-joints can be virtually tested under the same loading case by this modelling method to understand the effects of weave variations on the mechanical performance.

Table 1
Elastic properties of constituent materials (Units: GPa for moduli)

Table 2
Strengths of constituent materials (Units: MPa) currently available), as the responses are not very sensitive to them if they are in reasonable Load-displacement responses for the conformal (left) and voxel (right) model with different interface stiffness.Conformal: from 10 4 MPa/mm to 10 9 MPa/mm; response converged at k=105MPa/mm; Voxel: from 10 3 MPa/mm to 10 6 MPa/mm Table 3 Interface properties used in the 3D woven T-joint cohesive models