## Abstract

This paper is concerned with predicting the progressive damage and failure of multi-layered hybrid textile composites subjected to uniaxial tensile loading, using a novel two-scale computational mechanics framework. These composites include three-dimensional woven textile composites (3DWTCs) with glass, carbon and Kevlar fibre tows. Progressive damage and failure of 3DWTCs at different length scales are captured in the present model by using a macroscale finite-element (FE) analysis at the representative unit cell (RUC) level, while a closed-form micromechanics analysis is implemented simultaneously at the subscale level using material properties of the constituents (fibre and matrix) as input. The *N*-layers concentric cylinder (NCYL) model (Zhang and Waas 2014 *Acta Mech.* **225**, 1391–1417; Patel *et al.* submitted *Acta Mech.*) to compute local stress, srain and displacement fields in the fibre and matrix is used at the subscale. The 2-CYL fibre–matrix concentric cylinder model is extended to fibre and (*N*−1) matrix layers, keeping the volume fraction constant, and hence is called the NCYL model where the matrix damage can be captured locally within each discrete layer of the matrix volume. The influence of matrix microdamage at the subscale causes progressive degradation of fibre tow stiffness and matrix stiffness at the macroscale. The global RUC stiffness matrix remains positive definite, until the strain softening response resulting from different failure modes (such as fibre tow breakage, tow splitting in the transverse direction due to matrix cracking inside tow and surrounding matrix tensile failure outside of fibre tows) are initiated. At this stage, the macroscopic post-peak softening response is modelled using the mesh objective smeared crack approach (Rots *et al.* 1985 *HERON* **30**, 1–48; Heinrich and Waas 2012 *53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Honolulu, HI, 23–26 April 2012*. AIAA 2012-1537). Manufacturing-induced geometric imperfections are included in the simulation, where the FE mesh of the unit cell is generated directly from micro-computed tomography (MCT) real data using a code Simpleware. Results from multi-scale analysis for both an idealized perfect geometry and one that includes geometric imperfections are compared with experimental results (Pankow *et al.* 2012 *53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Honolulu, HI, 23–26 April 2012*. AIAA 2012-1572).

This article is part of the themed issue ‘Multiscale modelling of the structural integrity of composite materials’.

## 1. Introduction

Three-dimensional textile composites (3DTCs) are a relatively new class of ‘out-of-autoclave’ technology products that have the potential to replace conventional pre-impregnated (prepreg) tape laminates in terms of their excellent mechanical properties but at significantly reduced manufacturing cost and time. 3DTCs display high damage tolerance, high impact resistance and a distinct cost advantage and hence will probably see widespread use in the aerospace, automobile and defence industries due to their better structural properties tailoring capability. These composites can be created from a three-dimensional (3D) weaving or braiding process by arranging the fibre tows into complex dry preform structures, following which a resin is applied to the 3D preform to create the composite material. A significant amount of work has been done on the automated robotic approach for processing of textile preforms [1]. Two distinct types of 3D woven textile composites (3DWTCs), a layer-to-layer angle interlock 3DWTC and a Z-fibre orthogonal interlock hybrid 3DWTC, are popular and both have a great advantage over laminated materials. Angle-interlock 3D woven structures are common in order to create much thicker woven preforms. In the interlock structures, yarns can be woven from one layer to another and then back to the original layer to lock adjacent layers to each other. In complex interlock structures, yarns may be woven at specified points into several layers in order to join multiple layers. In this paper, Z-fibre orthogonal interlock hybrid 3DWTCs are studied in detail and analysed for progressive failure analysis.

Three-dimensional, orthogonal woven composites are a class of 3DWTC composite structures containing a set of fibre tows spreading in all three mutually perpendicular directions (*x*-, *y*- and *z*-axes). These structures display high resistance to layer delamination, as the latter is a common problem in conventional laminated composites. 3DWTCs are manufactured by laying up the warp and weft fibre layers and consequently running a Z-fibre in the thickness direction to bind the in-plane layers together. A matrix material is impregnated into these complex woven lay-ups and cured under specified conditions to obtain a woven solid structure. The weaving process followed ensures high consistency of the preform architecture and minimal waviness of the in-plane fibres in both the warp and weft directions. Several experimental studies [2–8] have shown that there is very low waviness of the warp and weft directional tows, which allows enhanced mechanical properties to be achieved. This category of composites is characterized with practically straight and well-aligned in-plane fibres. As a consequence, this type of textile composite shows significant improvement in in-plane elastic moduli and strengths than composites reinforced with more conventional 3D angle-interlock weaves usually produced by traditional two-dimensional (2D) weaving machines.

The hybridization of the composite material can be tailored for the specific needs of the structure by varying the fibre type and configuration of the weaving layers. Although there are many forms of 3D orthogonal woven composites with the variations of Z-fibre depth in the thickness direction, we will focus here mainly on one configuration, which is described in the following section. This hybrid architecture is provided by Textile Engineering and Manufacturing (T.E.A.M.), Inc., Woonsocket, RI, USA.

## 2. Material system

This paper discusses a ‘hybrid’ Z-fibre orthogonal interlock textile architecture 3D woven composite (H3DWC), where ‘hybrid’ refers to different constituent fibres, including IM-7 carbon, S-2 glass and Kevlar that are integrally woven into a single preform. The carbon fibres are used because of their high specific strength and stiffness, the glass for its relatively low cost (high strength per unit cost) and Kevlar for its high resistance (characterized by its toughness) to failure. Material properties for each fibre type and the matrix are provided in table 1.

## 3. Hybrid 2.5D woven composite configuration

A series of warp and weft fibres run in-plane throughout the panel with little or no undulation. A set of Z-fibres run in the direction of warp fibres and are drawn from bottom to top to bind all the layers together. Z-fibres are usually inserted between the spaces of the warp fibres. This architecture contains carbon layers for both outer surfaces (four layers of carbon on each side) and nine layers of glass in the middle. This hybrid configuration is a symmetric architecture that consists of orthogonally oriented in-plane carbon, and the glass tows for a total of 17 layers and referred to as ‘thick symmetric’. The Z-fibres run half the thickness in the two layers that bind the layers together and feature a distinct architecture compared with the through thickness features in other 3D architectures. These types of woven composites are categorized as 2.5D thick symmetric configuration for Z-fibres running half the thickness. The details of the configuration are shown schematically in figure 1, which provides some details of the Z-fibre path.

The planes that contain Z-fibres alternate along the weft direction with a period of two to achieve a repeating textile pattern. This configuration is infused with SC-15 polymer matrix material using a vacuum-assisted resin transfer moulding (VARTM) process. The geometric characteristics, unit cell dimensions and layer breakdown for this configuration are summarized in table 2.

The percentage content of the carbon and glass in this architecture is calculated by dividing the number of carbon and glass layers by the total number of layers of the composite. The image of the polished surface as shown in figure 2 demonstrates high straightness of the in-plane tows and high uniformity of the geometry. For weft tows, the waviness is highest in the first, fifth, sixth and ninth layers from the top, as caused by the Z-fibres, and is clearly visible in the figure. The fibre tow waviness and the crookedness are further examined with optical microscopy to characterize the textile architecture and geometric imperfections. The microscopy studies are discussed in §5.

## 4. Fabrications and manufacturing-induced imperfections

The hybrid orthogonal Z-fibre textile preforms are infused with SC-15 epoxy resin using a VARTM process to form a solid panel. VARTM is adapted from traditional resin transfer moulding (RTM) by replacing the upper half mould with a vacuum bag to enhance the impregnation of the fibre reinforcements. Details of the VARTM technique and fabrication process are provided in [9–12]. VARTM offers distinct advantages over RTM, including lower tooling cost, shorter mould time and ability to manufacture large structural components.

### (a) Geometry imperfections in the hybrid 2.5D textile composites

The Z-fibre orthogonal interlock hybrid 2.5DWCs do experience geometric distortion of the textile architecture during the fabrication process due to their weaving practice along the thickness direction. In this architecture, the Kevlar Z-fibre tows are used as weaves that are woven through multiple layers, showing small undulations along the weaving path. Compared with the rigid warp and weft tows running straight, the Z-fibre shows much more compliance, and, therefore, it is easily affected by the fabrication process, for example the tension exerted on the fibre tows during the weaving process and the mould pressure applied during the curing process. In the VARTM process with a single-sided mould, atmospheric pressure is exerted on the textile preform through the vacuum bag covering, forcing the fibre tows to be settled in a new position that is different from the predesigned one. This manufacturing-induced geometric imperfection of the textile architecture in the as-fabricated composite panel is evident from the microcomputed tomography (MCT) analysis of the sample, as described in detail in §5.

It has been pointed out by Song *et al.* [13] that each manufacturing process is associated with a unique set of characteristics that result in a produced part deviating from the expected ideal geometry. The set of such deviations, which is unique to each manufacturing process, is termed the ‘manufacturing imperfection signature’. Obtaining the manufacturing imperfection signature of the textile composite is important to determine the damage characteristics, such as strength, strain to failure and fatigue life, which has been studied by many researchers [14–22]. The importance of incorporating the unintended geometric deviations of the woven fabric into a textile architecture-based finite-element (FE) model has been recently addressed by Zhang *et al.* [23]. A voxel-based numerical model is studied by [24] to predict the deformed shape of 3D woven fabric during weaving and to investigate the effect of geometrical imperfections.

## 5. Microcomputed tomography analysis of geometry imperfections

In order to obtain a thorough understanding of the microstructure of the cured composite, cross-sectional microscopic images are used to characterize the textile geometry. The 2D images at various cross sections along the length of the specimen are rendered and reconstructed into a 3D volume for the purpose of characterization, and the measured dimensions are used as inputs to the textile architecture-based FE models presented in §5e.

### (a) Microcomputed tomography scanning parameter settings

The high-energy MCT scanner Skyscan 1173 with a 40–130 kV X-ray source and a pixel size of 6–140 μm was chosen as it has a large sample chamber for scanning objects of 140 mm in diameter and 140 mm in height. Specimens were cut along the warp and weft directions into 25×50 mm blocks and polished using a multi-speed rotating grinder/polisher. For this sample, we could potentially run between 70 and 80 kV with 17 μm resolution. A schematic of different steps carried out during the scanning process and reconstructing the 3D images is shown in figure 3. The MCT image data are analysed in detail, as described in the following sections.

### (b) Determining the unit cell

Unit cell scanned image data were cropped down from the larger specimen volume data. Cross sections along the mid-plane of fibre tows were chosen selectively to capture the periodicity of the unit cell in both in-plane directions. The unit cell in the black box is indicated in both the warp and weft directions, as shown in figure 4*a* and summarized in table 2.

### (c) Determining fibre tow cross-sectional details

Fibre tow cross-section dimensions are measured from the unit cell image data using a measuring tool in both the warp and weft directions, as indicated in figure 4*b*. Both the warp and weft tows are assumed to be straight with rectangular cross sections and the measured fibre tow dimensions are summarized in table 3. Scanning electron microscopy images are taken on the cross section of each type of constituent fibre tow to investigate how the fibres are distributed inside the tow and to determine the average fibre volume fraction. The fibre volume fraction of each constituent tow is summarized in table 4.

### (d) Porosity analysis

A study is carried out by calculating the volume percentage of the pores in the unit cell by adjusting the grey-scale difference in the scanned 3D image data. The volume fraction of pores in the unit cell is calculated as 0.91851%, which is negligible. Hence, the unit cell is assumed to be filled with materials without any pores and the effect of the pores is not included in the damage modelling and failure analysis of the textile composites.

### (e) Summary of microstructure output details from microcomputed tomography analysis

Here, all the outputs from MCT analysis are summarized as follows: (i) measurement of unit cell dimensions, (ii) measurement of fibre tow cross-section details in both the weft and warp directions and (iii) porosity analysis. An idealized 3D computer-aided design (CAD) model of the thick symmetric unit cell woven textile composite is modelled using the commercial software package ABAQUS 6.12, as shown in figure 5, taking into account all the above input values without any geometric imperfections. We will refer to this model as model I in the remainder of this paper.

### (f) Modelling *in situ* geometry imperfections

The FE modelling of 3D woven textile composites is still a challenging task for the research community because of the complexity of the textile architecture of fibre tows and very thin layers of matrix between winding fibre bundles. The manufacturing-induced geometrical imperfections have a strong impact on the performance and mechanical behaviour of the structure in the physical loading scenario and, hence, must be included in FE modelling to accurately capture its effect on the prediction of progressive damage and failure of the composite. In this paper, an innovative approach of numerical modelling of imperfections is demonstrated where a mesh is generated directly from MCT data using a commercially available software package, Simpleware. The point to be noted here is that the generated mesh is the real replication of the *in situ* microstructure imperfections.

Simpleware deals with real scanned image data from MCT analysis. It can generate FE meshes from MCT data, which helps in measuring and modelling the real imperfections in the structure. Each fibre tow is segmented separately and assigned material properties with a well-defined fibre orientation. The thin layers of matrix between the undulated fibre tows are meshed explicitly and assigned the material properties of SC-15 epoxy. In this novel and unique approach, modelling and meshing the *in situ* geometrical imperfections are possible in a detailed and efficient manner. These real measured imperfections can be taken as a reference to introduce imperfections to an idealized perfect geometry. The identified unit cell of 2.5D thick symmetric configuration is meshed explicitly using the features of Simpleware, as shown in figure 6. We will refer to this model as model II in the remainder of this paper.

A unidirectionally aligned fibre-reinforced composite is the most basic and simple microstructure for an aerospace composite material. It can be used to represent a fibre tow within a textile composite. Assuming that the fibres are infinitely long and the plane transverse to the fibre direction is statistically isotropic, the effective response is transversely isotropic for a fibre tow. The matrix material is an isotropic elastic damaging solid, and its nonlinear response is modelled using a modified *J*_{2} deformation theory of plasticity through a secant modulus approach, as discussed in §6. As a result, the composite effective stress versus strain response is extended to the nonlinear regime by substituting secant moduli into the stiffness tensor of the material system.

## 6. Matrix nonlinear constitutive model

Polymer matrix materials exhibit a nonlinear stress versus strain response prior to failure and when within a composite. Such nonlinear behaviour is attributed to matrix microdamage due to microcracking [25–28]. The nonlinear matrix response is modelled using a modified *J*_{2} deformation theory of plasticity, as formulated in §6a. The evolution of nonlinearity accounts for progressive deterioration of the material stiffness, however the tangent stiffness tensor still remains positive definite. The coalescence of matrix microdamage finally results in matrix macroscopic cracking, followed by a post-peak strain softening regime. In this instance, the positive definiteness of the matrix stiffness tensor is lost, which is categorized as a failure. The post-peak strain softening behaviour is modelled through the smeared crack approach (SCA) to capture matrix macroscopic cracking, as described in §6b.

### (a) Modelling microdamage in a polymer matrix

It has been shown by Sicking [29] and supported by experimental data from Lamborn and Schapery [30] that a polymer matrix exhibits limited loading path-independent behaviour, surmised through combined tension–torsion tests. Hence, for such a material, the state of stress can be uniquely determined from the state of strain through a secant modulus. It is further assumed that the evolution of damage is an irreversible process; therefore, once the matrix stiffness tensor is degraded due to microdamage, it cannot be recovered. Such behaviour suggests that a modified *J*_{2} deformation theory of plasticity can be employed to model the nonlinear stress–strain response. In the present constitutive model, the degrading secant moduli are used to compute the material stiffness tensor during unloading also.

In order to use a uniaxial stress–strain response to determine the material responses under multi-axial loading, two equivalent variables, the equivalent stress, *σ*^{m}_{eq}, and the equivalent strain, *ϵ*^{m}_{eq}, are introduced and related through a secant Young modulus, *E*^{m}_{s}, as . The two equivalent quantities are related to the stress and strain components as [31]
6.1and
6.2where *ν*^{m}_{s} is the matrix secant Poisson ratio defined by
6.3and *E*^{m}_{e} and *ν*^{m}_{e} are the matrix elastic Young modulus and elastic Poisson ratio, respectively. In this paper, the matrix nonlinear stress–strain relation is characterized using an exponential relation as
6.4where is the proportional limit determined from a uniaxial tensile test and *K*_{1} and *K*_{2} are the two material constants to fit the experimental nonlinear stress–strain response. As discussed in [28], the matrix *in situ* response can be characterized through a tensile test on a ±45° symmetric laminate. Different nonlinear responses have been reported for the matrix within a carbon tow from that within a glass tow [32]. Table 5 summarizes the matrix nonlinear properties used in the two-scale micromechanics model (discussed in §7) for computing the response of each constituent tow. In addition, the pre-peak nonlinear properties of the surrounding matrix in the macroscale model is assumed to be the same as that used for the glass tow.

### (b) Modelling failure progression in a polymer matrix using the smeared crack approach

The accumulation of matrix microdamage leads to the initiation of matrix macroscopic cracking, followed by a post-peak strain softening response. In this study, the evolution of matrix failure was modelled using the SCA, which was originally developed by Rots *et al.* [33] to model crack propagation and fracture in concrete. In the SCA, it is assumed that distributed cracks are ‘smeared’ out over a certain width within the FE such that the effect of progressive cracking is represented by macroscopic strain softening in a continuum scheme, as shown in figure 7. To restore mesh objectivity, a characteristic length is introduced such that the total amount of energy dissipated during failure in a continuum element is equal to the fracture toughness defined for a cohesive element of the same size. The fracture toughness, or the critical energy release rate, *G*_{C}, is defined by the area under the traction-separation law that dictates the cohesive behaviour of crack propagation (figure 8*b*) as
6.5where *u* is the sum of crack displacements within the fracture zone, as schematically shown in figure 7*a*. In the SCA, *u* represents the crack strain acting across a certain width within a finite element, denoted as the crack band width, *h*, which is shown in figure 7*b*. Assuming that all the cracks are uniformly distributed over the crack band, and since *u* is the accumulation of all the crack strains over the fracture zone, it follows that
6.6If *g*_{c} represents the area under the softening branch of the stress–strain response (figure 8*a*), then substituting equation (6.6) into equation (6.5) results in
6.7

Therefore, the strain-based description for a softening material is related to the displacement-based traction-separation laws through the characteristic length, *h*. In an FE setting, *h* is chosen based upon the element type, element size, element shape and the integration scheme [33]. Typically, the length of the element projected onto the crack normal is used as a characteristic element length, as shown in figure 9.

#### (i) Formulation

The formulation of the SCA for matrix cracking is similar to the one presented in [33,34]. In the pre-peak regime, the material is governed by a standard continuum theory such as elasticity or plasticity. When failure is initiated, it is assumed that the total strains, ** ϵ**, are decomposed into continuum strains,

*ϵ*^{co}, and cracked strains,

*ϵ*^{cr}, as 6.8

In a continuum scheme, *ϵ*^{co} can be further decomposed into elastic, plastic and thermal strains, if needed. In this study, the matrix nonlinear response is modelled as a degrading secant solid; hence, the secant strains at the onset of failure initiation are used for *ϵ*^{co}. Here, ** ϵ**,

*ϵ*^{co}and

*ϵ*^{cr}are presented in global coordinates.

Figure 10 shows the crack morphology in three dimensions. At the crack interface, there exist three relative displacements between the crack faces: one is the crack opening displacement, and the other two are the crack sliding displacements. The subscripts *n* and *t* are used to designate directions normal to the crack and tangential to the crack, respectively. The key to the SCA is to embed cracks into a continuum; hence, the crack opening displacement is represented by a local normal crack strain, , and the two crack sliding displacements are replaced by two local shear crack strains, and . These local crack strains are transferred to the global coordinates through
6.9where *e*^{cr} is a vector that contains local crack strains and ** T** is a 6×3 transformation matrix depending on the crack orientation (appendix A). Similarly, the interface stresses at the crack interface,

*s*^{cr}, can be transferred to the global stress state,

**, through 6.10**

*σ*The crack interface stresses are related to the local crack strains through
6.11where *D*^{cr} is the crack interface stiffness matrix that dictates the failure evolution in the post-peak strain softening regime. For a single crack in a 3D solid, *D*^{cr} can be expanded as
6.12where *D*_{c} is the secant stiffness across the crack interface due to crack opening and *G*_{c1} and *G*_{c2} are the two secant shear stiffnesses governed by crack sliding. These quantities identify the modes of failure and are related to the traction-separation laws with a characteristic length scale. The off-diagonal terms are assumed to be zero, indicating that there is no coupling between the normal and shear crack components. The coupling between the crack shear and opening displacements, known as crack dilatancy, has been extensively studied by Bažant & Gambarova [35], Walraven [36], Walraven & Reinhardt [37] and Gambarova & Karakoç [38].

The constitutive relation for a continuum is
6.13where *D*^{co} is the continuum stiffness tensor. In the current study, *D*^{co} is computed using the secant Young modulus, *E*^{m}_{s}, and the secant Poisson ratio, *ν*^{m}_{s} (see equation (6.3)) at the onset of failure initiation. Combining equations (6.8)–(6.10) and (6.13) results in
6.14Finally, the relation between the total stress and total strain in the post-peak regime is computed as
6.15

Details of the numerical implementation of the SCA within an FE model including the iteration scheme to determine the crack strains are discussed in [31,34].

#### (ii) Determination of the crack constitutive relations

In this study, a one-dimensional (1D) uncoupled traction-separation law is employed; consequently, equation (6.12) becomes
6.16where *E*^{cr} is the secant crack modulus resulting from normal crack strain (mode I type of failure), while and are associated with shear crack strains (mode II type of failure). In this study, secant crack stiffness is used so that the softening response follows the traction-separation law exactly, as shown in figure 11. To prevent crack healing from happening, it is required that
6.17Hence, once the crack stiffness is degraded, it cannot be recovered. The loading and unloading responses during the evolution of the failure process are specified in figure 11.

In the experiment, matrix cracking is observed in regions of predominant tensile stress, as presented in [39]. In a monolithic material, cracks are likely to grow under pure mode I conditions as this mode of failure is energetically favourable. In H2.5DTCs, the cracks in the surrounding matrix may be subjected to mixed-mode loading due to the presence of fibre tows. However, in this study, it is assumed that the tensile cracks grow under pure mode I conditions, oriented with the maximum principal stress plane. For a given stress state, the principal stresses, *σ*_{1}, *σ*_{2} and *σ*_{3}, and the corresponding principal axes, *n*_{1}, *n*_{2} and *n*_{3}, are first computed. The maximum tensile stress criterion is used for tensile crack initiation as
6.18where *σ*^{cr} is the cohesive strength, or the critical stress, in the mode I traction-separation law. In practical applications, it is further assumed that, once the crack is initiated, the crack orientation, determined from the principal stress directions is fixed during the failure evolution. It is further assumed that the crack interface is traction free in both normal and shear directions during and after fracture energy dissipation. The two crack shear moduli, *G*^{cr}_{1} and , are degraded as a function of , indicating that the cracks grow under mode I dominated conditions. It is possible that the crack shear moduli are degraded with respect to or , and a mixed-mode traction-separation law could be introduced to ensure that the shear failure evolves under mode II conditions [40]. However, such a complicated failure mechanism requires further study of cracks progressing at the microscale, and this aspect is not considered in this study. The fracture properties of the surrounding matrix used in the macroscale model are given in table 5.

## 7. Fibre tow constitutive model

A curved fibre tow within a textile composite can be represented using a multi-layered concentric cylinder model. The fibre packing inside the tow can be assumed as (i) a hexagonal array in the transverse plane where the fibres are of identical cross section or (ii) a random array where the fibres are of different diameters but the volume fraction is preserved constant. In both cases, the fibres are assumed to be infinitely long and the effective response of the composite is macroscopically homogeneous and transversely isotropic, requiring five independent constants to form the composite stiffness tensor. Though the choice of these elastic constants is not unique, the axial modulus, , the axial Poisson ratio, , the axial shear modulus, , the plane-strain bulk modulus, , and the transverse shear modulus, , are used throughout the paper. Therefore, the stiffness tensor for a transversely isotropic material can be written in terms of these elastic constants as 7.1

Other important constants, including the transverse modulus, , and transverse Poisson ratio, , can be computed as 7.2and 7.3where

In this study, the fibre is assumed to be linear elastic, transversely isotropic, with ‘1’ designating its longitudinal direction. Its stiffness tensor, *C*^{f}, can be written in terms of fibre elastic properties as equation (7.1) by replacing the superscript ‘c’ with ‘f’. The fibre tows exhibit nonlinear stress–strain responses due to the evolution of matrix microdamage (discussed in §6a) developing at the microscale. Such nonlinear responses can be captured using the two-scale micromechanics model recently developed by Zhang & Waas [41]. This two-scale micromechanics model is extended to any number of matrix layers, in general, by Patel *et al.* [42], and is henceforth called the *N*-layer concentric cylinder (NCYL) model. The implementation of this NCYL model for a fibre tow in a multi-scale computational framework for H2.5DTCs is presented in §7a. In the present NCYL model, the inner fibre core is surrounded by (*N*−1) matrix layers. During the process of damage evolution, the damage occurring in the discrete matrix layers is identified based on the calculation of the strain and stress fields, whereas the undamaged matrix layers retain their virgin properties. It is important to emphasize the fact here that, in the 2CYL model, the entire volume of matrix material was degraded by the same amount, as it is not possible to restrict damage to a particular region of matrix volume. As a result, the model overpredicts the stiffness and the nonlinear stress–strain response. This limitation is eliminated by the NCYL model and the radius of the discrete damage zone can be controlled as an input variable to this model. The homogenized elastic properties of the fibre tow (, , , and ) are functions of the degraded properties of the damaged matrix layers and the virgin properties of the undamaged matrix layers. This is very crucial for progressive damage analysis, which is based fully on an analytical solution.

In the experiment, catastrophic failure modes are reported, including fibre tow breakage, tow splitting and transverse cracking, resulting in a loss of load-carrying capability at the macroscale, followed by post-peak strain softening in the stress versus strain response [39]. Since the positive definiteness of the tangent stiffness matrix is lost, the aforementioned two-scale model will provide mesh-dependent results in an FE framework if no characteristic length is introduced. Therefore, the SCA established in §6b is suitable for modelling the tow failure behaviour, as presented in §7b.

### (a) Progressive damage associated with the pre-peak nonlinear response

#### (i) Micromechanics (NCYL model)

The micromechanics analysis at the fibre and matrix scale is able to predict the fibre tow effective properties using the NCYL model. In this model, the composite material is represented by an inner fibre core and an outer matrix annulus. The matrix annulus is divided into (*N*−1) layers of concentric cylinders, where the ratios of concentric matrix cylinder radii are kept constant in general. It is assumed that the size of the concentric pair varies such that the entire volume of the composite can be filled with the cylinders, while the ratio of the fibre radius to the matrix outer radius is kept constant to maintain the fibre volume fraction. A representative *N*-layer concentric cylinder unit with fibre radius *a* and matrix outer radius *b* is shown in figure 12, and the resulting fibre volume fraction is *V*_{f}=*a*^{2}/*b*^{2}

The two-layer fibre–matrix analytical model [41] has been extended to *N*-layers concentric cylinders [42] to analyse the stress and strain fields for all constituent matrix layers. In order to resolve the local fields in the constituent materials at the micromechanics scale, a composite cylinder with a multiple number of matrix layers is considered in such a way that each layer corresponds to a different volume. The detailed description of the NCYL micromechanics model and its implementation to multi-scale analysis for progressive damage and failure analysis of the composite structure is explained in [42].

The key to the micromechanics model proposed by Zhang & Waas [41] is to relate the fibre tow strains (the strain applied on the outer matrix cylinder) to the local matrix strain fields through a 6×6 transformation matrix, ** F**, as,
7.4The

*F*

_{ij}components can be computed by imposing a single non-zero composite strain on the fibre–matrix cylinder and solving the resulting matrix strain fields. In particular, the axial properties, including the axial tensile ( and ) and axial shear (), are computed through a two-phase concentric cylinder model (CCM) [43], which is subsequently used for the computation of

*F*

_{i1},

*F*

_{i4}and

*F*

_{i5}. The rest of the components in the

**matrix are determined via an extended three-phase generalized self-consistent model (GSCM) [44], which also gives the composite transverse properties, and . The computations of tow effective properties and each component in the**

*F***matrix for**

*F**N*=2,3 and 4 cylinders are provided in [42]. It should be noted that the proposed micromechanics model is based upon the homogenization techniques that were originally used to compute the composite effective moduli, hence both the fibre and matrix are assumed to be linear elastic. When the matrix stiffness is reduced due to microdamage, the nonlinear response of the matrix is modelled through a secant modulus approach, in which the matrix elastic properties are replaced with the corresponding secant moduli.

#### (ii) Multi-scale computational framework for nonlinear fibre tow

In the two-scale modelling methodology established in this study, the macroscale, fibre tow level analysis is conducted by using effective homogenized properties to compute the local stress and strain fields in the fibre tow. Simultaneously, it is intended to carry out the subscale analysis, at the fibre and matrix level, using the micromechanics model presented in §7a(i), in which the constituent stress and strain fields are calculated in closed form.

The commercially available finite-element software ABAQUS (v. 6.12) is used for the macroscale FE analysis, and the NCYL micromechanics model at the subscale is implemented at each integration point of the macroscale, using a user-defined material subroutine, VUMAT (for the ABAQUS/Explicit solver). This subroutine is called at each integration point at each increment, and the material constitutive law is updated through user-defined options. At the start of each increment, the material states, i.e. the stress–strain and solution-dependent state variables from the previous equilibrium step and the strain increments in the current step, are passed into the VUMAT through the ABAQUS solver.

#### (iii) Implementation of the two-scale model

The strains at each integration point in the FE analysis of the homogenized fibre tow are applied to the subscale micromechanics NCYL model. These integration point strains can be treated as the effective composite strains that are applied on a discrete fibre–matrix microstructure. The constituent strain fields therefore can be computed in closed form by knowing the globally applied strains through the micromechanical analysis discussed in §7a(ii).

However, it should be noted that the resulting matrix strain fields vary in space; hence, the matrix equivalent strain, computed using equation (6.2), has spatial variation as well. In the current fully analytical computational scheme, it is hypothesized that the composite nonlinear behaviour can be characterized using two scalar variables that are related to the matrix equivalent strain. This idea is similar to the mean-field theories in which the average values of the strain fields are used to determine the matrix nonlinear progression. In this study, the two scalar variables that govern the evolution of matrix microdamage are defined based upon the maximum and average value of the square of the matrix strains at the fibre–matrix and matrix–matrix interfaces, respectively, as
7.5and
7.6where is the matrix strain at the fibre–matrix and other matrix–matrix layer interfaces. Physically, the average term is dominant in the matrix strain field when the fibre volume fraction is low, while the maximum value dominates the result for the high fibre volume fraction. Therefore, a weight function of *V*_{max} and *V*_{avg} can be written as
7.7where *n* is dependent on the fibre-to-matrix stiffness ratio such that the effect of constituent properties can be accounted for. In this paper, it is assumed that
7.8Consequently, two matrix equivalent strains can be computed: one is based upon the weight function in equation (7.7), while the other is based upon the average value in equation (7.6) as
7.9

Once the matrix equivalent strain is resolved, the matrix stiffness tensor is degraded as a secant solid according to the nonlinear damage model presented in §6. Consequently, if matrix microdamage occurs, the stiffness of the subscale microstructure is reduced based upon the proposed secant modulus approach. The subscale stiffness tensors are subsequently used to update the global stiffness and stresses in the macroscale FE analysis. Therefore, the influence of the matrix microdamage at the subscale manifests as the degradation of the global stiffness of the composite. In order to verify the proposed method to compute the composite nonlinear response, a 3D FE model is used to assess the accuracy of the predictions as explained in [42]. A flowchart of the multi-scale approach followed in this study is shown in figure 13.

### (b) Tow failure mechanisms: post-peak strain softening response

In H2.5DWCs, due to the heterogeneity of the microstructure and the complexity of the local fields, fibre tows exhibit multiple failure modes, including tensile tow breaking and tow splitting in transverse and shear responses. In some instances, the fibre tow can be delaminated from the surrounding matrix at high strain rates, as shown by Pankow *et al.* [45] through a split Hopkinson pressure bar test. Shear bands are observed when the textile composite is subjected to the through-the-thickness compression [10,45]. Cracks that grow along the transverse normal, axial shear and transverse shear directions are considered as matrix failure modes, dominated by the strength and toughness of the polymer matrix material.

Since the fibre tow pre-peak nonlinear response is attributed to matrix microdamage, no macroscopic damage criterion is required to drive the nonlinear damage progression. However, multiple catastrophic failure modes are observed in the experiment, including tensile failure due to fibre breakage and fibre tow splitting. These modes of failure result in a loss of load-carrying capability of composite structure, followed by a post-peak strain softening response. Since the positive definiteness of the material tangent stiffness matrix is lost in the softening regime, the FE analysis will provide mesh-dependent results if no characteristic length is introduced. As a result, the aforementioned two-scale model has to be supplemented by a suitable mesh objective approach for modelling the post-peak softening response.

In this study, fibre tow breaking and splitting are considered as the two major failure modes for H2.5DWCs subjected to tensile loading, which is evident from the experimental investigation discussed in [39]. In addition, when the tow breakage is initiated in axial tows due to tensile stress, the crack plane is assumed to be aligned perpendicular to the fibre direction. On the other hand, the crack plane is assumed to be parallel to the fibre direction when the tow splitting criterion is reached due to transverse tensile loading. Details of the implementation of the SCA within a FE framework are demonstrated in §7b(i).

#### (i) Implementation of the smeared crack approach to model the failure response of a single fibre tow

The SCA formulated in §6b is used to model the failure progression of a single fibre tow, which includes tensile breakage in the axial direction and tow splitting in the transverse direction. It is further assumed that, when the critical stress is reached, the crack plane is aligned perpendicular to the fibre direction for tow breaking and the crack orientation transformation matrix, ** T**, is determined as given in appendix A. Tow splitting occurs due to matrix cracking inside the fibre tow subjected to transverse and shear loading. In this work, it is assumed that, when the transverse critical stress is reached, the crack plane is aligned parallel to the fibre direction for tow splitting. Therefore, the crack orientation transformation matrix,

**, is determined as given in appendix B. When two-piece failure occurs in tensile loading, the recorded load drops significantly during the experiment, indicating a considerable amount of fracture energy dissipation.**

*T*In this study, it is assumed that, when the tensile failure mode occurs, the fracture energy is dissipated completely. The fibre tow failure properties, including the critical stress and fracture toughness for each mode of failure, used in the macroscale FE model are given in table 6. The tensile strength values for carbon and glass tows are taken from [46,47], respectively. It is worth mentioning that characterizing the failure progression of fibre tows within H2.5DWCs is critical to understand the progressive failure response of this class of materials. Aspects including the analysis of tow splitting due to shear loading, shear band formation and subsequent energy dissipation during failure are relegated to a separate study.

## 8. Results and discussion

Stress–starin responses, linear elastic moduli, Poisson ratios, peak strength and failure strains have been experimentally determined in [39]. The tensile test simulations are carried out for model I (idealized perfect model with no imperfections) and model II (imperfect model with *in situ* geometric imperfections generated by Simpleware) simultaneously in both the weft and warp directions and the simulation results are compared with the experimental results as summarized in the following sections. Model I had nearly 50 000 elements and took 20 h to run on a high-performance computing system using 32 CPUs, whereas model II had nearly 590 000 elements and took 2 days to run on the same computing system using 32 CPUs.

### (a) Elastic and strength properties

Considerable gradual linear growths of the elastic moduli under loading in the weft and warp directions are observed for strains from 0.1% to 0.5% (figure 14). The linear growth is followed by a slow gradual reduction of the moduli with further increase of strains. The nonlinear behaviour is attributed to a combination of matrix microdamage and fibre splitting. The latter appears to be a minor cause contributing less to the modulus versus strain variations. The stress–strain response, as shown in figure 14*a*, exhibits a higher degree of linearity in the weft direction due to the fact that loading occurs directly on the fibre tows situated in the outermost surface. The warp direction shows more of a progression in failure of the material as shown in figure 14*b* and deviations from linearity occur relatively early in the loading regime. However, due to the fact that different fibre tows will have different stresses at the same externally applied displacement indicates that the carbon may fail earlier than the glass and the failure envelope is progressive due to the hybrid materials before two-piece failure. Numerical predictions using the two-scale multi-scale method showed excellent agreement with experimental data for in-plane elastic moduli in both the weft and warp directions for the idealized perfect geometry without any imperfections (model I). After including the *in situ* geometrical imperfections and the fibre crookedness in the FE model (model II), the predicted elastic moduli reduced the respective model I values (e.g. warp-directional modulus 34.12 GPa and weft-directional modulus 39.95 GPa) by only 2.2% in the warp and 2.3% in the weft directions. As we can see, there is a reduction in stiffness after including the geometric imperfections. For this hybrid configuration, the weft direction is stiffer than the warp direction due to the fact that there is one additional glass fibre tow layer in the weft direction, as summarized in table 7.

The weft direction tensile strength was found to be 467.61 MPa and the warp direction strength was 446.81 MPa, as reported in table 8 for model I; the former is 4.7% larger due to a 5.8% higher fibre weight fraction in the weft direction. In model II, because of the imperfections of fibre alignment in both the in-plane and thickness directions, there is significant fibre tow splitting observed in the transverse and Z-fibre tows. The predicted tensile strength reduced the respective model I values (e.g. warp-directional strength 403 MPa and weft-directional strength 449.7 MPa) by only 9.8% in the warp and 3.8% in the weft directions. As observed, there is a degradation of tensile strength after including the geometric imperfections, and predictions are closer to the experimental values. The global response tends to be highly nonlinear and matches well with experimental results for model II as shown in figure 14*b*, which demonstrates the effect of including manufacturing-induced imperfections in the progressive damage and failure analysis of these textile composites.

### (b) Progressive damage during quasi-static tensile loading

The multi-scale methodologies described in §§6 and 7, which combine modelling pre-peak nonlinearity using the NCYL secant stiffness method and the post-peak strain softening response using SCA, are used together to conduct a thorough investigation of the damage and failure mechanisms in the hybrid textile composite under consideration. The progressive damage and failure responses of textile representative unit cells (RUCs) in both the weft and warp directions subjected to uniaxial quasi-static tensile tests are shown in figures 15–18 for imperfect unit cells with geometrical imperfections (model II). In these figures, the progressive failure statuses at different percentages of the failure strength are shown as contour plots of progressive failure flags of the elements.

First, a delayed damage initiation has been observed in the weft direction loading. First, matrix cracking in the warp direction loading is observed at 0.45% strain, which is much earlier than weft loading (observed at 0.67% strain). The damage initiation for the loading in the weft or warp directions occurs in the range of applied strain 0.4–0.7% (figures 15 and 17, showing the data for weft and warp direction loading). This range of damage initiation strains is significantly higher than the respective ranges observed for other carbon/epoxy textile composites and is at the lower end of the typical damage initiation strain range (0.6–0.7%) for conventional cross-ply prepreg tape laminated composites.

The damage in these composites starts near the location of the Z-fibre tows and at the edges of the in-plane fibre tows oriented transversely to the loading direction. Because of this, there is a significant amount of local disbonds and matrix cracking occurring in the case of warp direction loading, as shown in figure 18. This failure event leads to the development of transverse cracks inside the in-plane and vertical tows (referred to as ‘tow splitting’) during the progressive loading process; these are followed by fibre tow breakage in the axial direction at a very late stage of loading and finally cause the ‘two-piece’ failure. A comparison study is being conducted from the simulation results for progressive damage and failure, for both weft and warp direction loading. The study is based on figures 15 and 17 and is summarized in table 9.

Damage initiation and evolution in both the weft and warp directions are investigated numerically using the NCYL multi-scale framework and the matrix microdamage model based on *J*_{2}-deformation theory of plasticity. Previous experimental studies [48–50] discuss the damage progression in 3D woven composites in a detailed manner, which correlates well with the sequence of damage events captured with our proposed model. Transverse crack kinetics within a homogenized fibre tow is explained in the literature [51] by bridging the micro-meso scale. We have implemented the physics of this failure mode in our smeared crack code and captured the transverse cracks. A micromechanics-based damage approach is applied to the macroscale unit cell model by Saleh *et al.* [50], which approximates the 3D woven composite as a cross-ply laminate without considering the effects of the binder yarns in the thickness direction and also neglects the void content and the geometrical imperfections induced due to the manufacturing process. Our modelling and analysis of unit cells includes *in situ* geometric imperfections and the predicted simulation results match closely with the sequence of damage events described by them. The scope of this work is to establish a computationally efficient progressive damage and failure analysis tool and the goal is achieved by using the NCYL micromechanics multi-scale model, which is based purely on an analytical solution. The recommended matrix damage model improvements will be the focus of a future study for off-axis tension and shear loading simulations.

Overall, the FE model results that include imperfections (model II) show a good prediction for the elastic stiffness in each case. Also, the model is able to capture the nonlinear stress–strain response to a close extent. It is noted that this model includes real *in situ* microstructure imperfections inside the unit cell and the multi-scale analysis is carried out to capture the progressive failure. The differences in global stress–strain responses between the idealized perfect unit cell (indicated as a blue line) and the one with real *in situ* microstructure imperfections (indicated as a black line) are visible in figure 14.

### (c) Uniaxial coupon-level tensile test for multiple representative unit cells

Test specimens of suitable dimensions are to be decided based on the number of RUCs, which should be sufficient to run the simulations and capture the main features of the progressive damage and failure characteristics. Nine RUCs (3×3) are assembled and studied for uniaxial coupon-level test simulations using the multi-scale analysis. A breakdown of the 3×3 RUCs for a thick symmetric (2.5D) woven textile composite is shown in figure 19 for illustration.

The coupon-level global stress–strain response is shown in figure 14*a* for the weft direction, which matches fairly well with the single RUC results. Also, the progressive damage and failure response of the coupon in the weft direction subjected to uniaxial tensile loading is shown in figure 20. In these figures, the progressive failure statuses at the failure strength are shown as the contour plots of progressive failure flags of the elements.

## 9. Summary and conclusion

In this study, hybrid 2.5D thick symmetric woven textile composites with symmetric fibre architecture at the top and bottom and two layers of Z-fibres running through the thickness direction are simulated for uniaxial tensile responses to determine the effect of hybridization and compute the effective stiffness and tensile strength for both the warp and weft directions, including progressive damage and failure. An MCT analysis is carried out to determine and measure real *in situ* microstructural imperfections. This MCT analysis helps to extract essential inputs (cross-sectional details of fibre tows in both the warp and weft directions and volume fraction of pores) for constructing RUCs. Based on these inputs, an idealized 3D CAD model of a unit cell is constructed and multi-scale analysis is carried out for progressive failure analysis. A subscale micromechanics 2CYL model is used to establish a multi-scale computational framework to predict the effective nonlinear response of a homogenized fibre tow. The influence of matrix microdamage at the microscale manifests as the degradation of the effective fibre tow stiffness at the macroscale through a secant modulus approach. Since fully analytical solutions are used for the subscale micromechanics analysis, the proposed method offers a lower computational cost and is suitable for large-scale progressive damage and failure analysis of composite structures. The multi-scale homogenization of the fibre tows is used as the building blocks for an explicit FE model of 3D woven textile composite RUCs. The linear elastic stiffness of this textile architecture matches well to the experimental results in both the warp and weft directions. Also, the predicted failure strength and the global nonlinear stress–strain responses are in good agreement with the experiments.

It is a challenging task to accurately model geometric imperfections, such as crookedness and misalignment of fibre tows, to study the effect of microstructure imperfections on the global response of a textile composite structure. An effective software tool, Simpleware, is used to generate a FE mesh from real MCT *in situ* data; this mesh is able to capture the microstructure details and provide a platform to carry out multi-scale analysis on a real RUC model. The progressive failure responses from both models, namely (i) an idealized unit cell with no imperfections (model I) and (ii) a Simpleware-generated model with *in situ* imperfections (model II), are compared to demonstrate the effect of including geometric imperfections on the overall behaviour of the composite. The coupon-level simulation is carried out to understand the size effect by a combination of multiple RUCs. Overall, the main objective of this paper is achieved by establishing a multi-scale method which is capable of predicting the effective linear elastic stiffness, the global nonlinear stress–strain response and the two-piece failure strength of hybrid 2.5D thick symmetric textile composites. The progressive analysis methodology can be extended to other configurations of woven and braided textile composites including the effect of hybridization and architectural variations. In this study, the multi-scale analysis combines the 2CYL micromechanics model for the pre-peak nonlinear response and the SCA to model the post-peak failure analysis. This combined computational approach is used in model I and model II. The approach in model II, which uses a mesh generated directly using MCT data, is novel and provides a means to capture unintended microstructural imperfections due to manufacturing process-related effects.

## Competing interests

We declare we have no competing interests.

## Funding

The authors are grateful for the financial support from the Army Research Laboratory, Aberdeen Proving Ground, MD, and the Army Research Office.

## Acknowledgements

The authors are grateful for the encouragement and support of Dr Asher Rubinstein, Program monitor at the U.S. Army Research Office, Research Triangle Park, North Carolina, USA.

## Appendix A. Calculation of the crack orientation transformation matrix for tow break and epoxy tensile failure

Assuming that the global coordinates are represented by the 1−2−3 triad while the local crack orientation is denoted by the 1′−2′−3′ triad, with the crack normal aligned with the 1′-direction, it follows that
A 1where
A 2and *a*_{ij}′*s* are the direction cosines governing the space vector transformation as
A 3

When the SCA is implemented for a fibre tow, the crack plane is assumed to be perpendicular to the fibre direction. Thus, the 1′-axis that defines the crack normal coincides with the 1-axis that denotes the fibre direction, and ** N** is reduced to
A 4

## Appendix B. Calculation of the crack orientation transformation matrix for the tow split

Assuming that the global coordinates are represented by the 1−2−3 triad while the local crack orientation is denoted by the 1′−2′−3′ triad, with the crack normal aligned with the 1′-direction, it follows that
B 1where
B 2and *a*_{ij}′*s* are the direction cosines governing the space vector transformation as
B 3

When the SCA is implemented for a fibre tow split due to transverse loading, the crack plane is assumed to be parallel to the fibre direction. Thus, the 1′-axis that defines the crack normal is perpendicular to the 1-axis that denotes the fibre direction, and ** N** is reduced to
B 4

## Footnotes

One contribution of 22 to a Theo Murphy meeting issue ‘Multiscale modelling of the structural integrity of composite materials’.

- Accepted March 18, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.