Multiscale Progressive Failure Analysis of 3D Woven Composites

Application of three-dimensional (3D) woven composites is growing as an alternative to the use of ply-based composite materials. However, the design, analysis, modeling, and optimization of these materials is more challenging due to their complex and inherently multiscale geometries. Herein, a multiscale modeling procedure, based on efficient, semi-analytical micromechanical theories rather than the traditional finite element approach, is presented and applied to a 3D woven carbon–epoxy composite. A crack-band progressive damage model was employed for the matrix constituent to capture the globally observed nonlinear response. Realistic microstructural dimensions and tow-fiber volume fractions were determined from detailed X-ray computed tomography (CT) and scanning electron microscopy data. Pre-existing binder-tow disbonds and weft-tow waviness, observed in X-ray CT scans of the composite, were also included in the model. The results were compared with experimental data for the in-plane tensile and shear behavior of the composite. The tensile predictions exhibited good correlations with the test data. While the model was able to capture the less brittle nature of the in-plane shear response, quantitative measures were underpredicted to some degree.


Introduction
Three-dimensional (3D) woven composites have emerged as attractive alternatives to traditional laminated composites by offering improved damage tolerance, customizability, and near net-shape manufacturing. This enabling technology allows structures with thick sections and/or complex geometries to be fabricated using a single preform that can then be infused with resin. Woven composites with new 3D architectures have predominantly been evaluated experimentally due to the significant challenges associated with modeling 3D woven geometries. Macromechanical modeling of these materials is particularly problematic because of their strong dependence on microstructure, coupled with their more general anisotropic behavior (compared to tape and 2D woven composites). As such, micromechanics and multiscale modeling are attractive for examination of these complex materials.

Geometric Characterization Based on X-ray CT and SEM
For the IM7/RTM6 3D woven composite considered in this investigation, X-ray images (13.4 µm resolution) were previously generated for individual pristine specim [30]. These data are from panel "1Z" in ref. [30], and readers are referred to that refer for more details regarding data generation and observed features. Note that the cou whose test data are presented herein came from this same panel. In this work, the X CT data were re-analyzed to identify and quantitatively measure various features o terest.
Based on features observed in the X-ray CT data, a number of key dimensions w identified and measured for use in the model RUC. Representative images are show Figure 2. The measured dimensions include the warp-tow width (WaW) and he (WaH) and the weft-tow width (WeW) and height (WeH). For sections of the binder aligned in the warp direction, the elliptical binder-tow shape is similar to the shap the warp-and weft-fiber tows. As the binder-tow transitions through the thicknes width narrows and its height increases. To approximate this behavior, two sets of bin tow dimensions were measured: the width and height in the flat warp-aligned reg (BWf and BHf, respectively) and near the mid-thickness (BWz and BHz, respectiv Statistical distributions of these eight dimensions were determined by repeatedly ext ing dimensions from different locations near the appropriate tow centerline and from ferent cross-sectional planes if possible. For example, the warp-tow width, WaW, ca found by measuring warp tows from images normal to the warp direction or norm the thickness direction. The mean values of the identified key dimensions were use the construction of the 3D woven RUC for the model described in Section 2.3.3. In future, automated techniques, such as those described in refs. [31][32][33][34][35], could be utilize speed up the geometry-characterization process.

Material Characterization 2.2.1. Geometric Characterization Based on X-ray CT and SEM
For the IM7/RTM6 3D woven composite considered in this investigation, X-ray CT images (13.4 µm resolution) were previously generated for individual pristine specimens [30]. These data are from panel "1Z" in ref. [30], and readers are referred to that reference for more details regarding data generation and observed features. Note that the coupons whose test data are presented herein came from this same panel. In this work, the X-ray CT data were re-analyzed to identify and quantitatively measure various features of interest.
Based on features observed in the X-ray CT data, a number of key dimensions were identified and measured for use in the model RUC. Representative images are shown in Figure 2. The measured dimensions include the warp-tow width (WaW) and height (WaH) and the weft-tow width (WeW) and height (WeH). For sections of the binder tow aligned in the warp direction, the elliptical binder-tow shape is similar to the shapes of the warp-and weft-fiber tows. As the binder-tow transitions through the thickness, its width narrows and its height increases. To approximate this behavior, two sets of bindertow dimensions were measured: the width and height in the flat warp-aligned regions (BWf and BHf, respectively) and near the mid-thickness (BWz and BHz, respectively). Statistical distributions of these eight dimensions were determined by repeatedly extracting dimensions from different locations near the appropriate tow centerline and from different cross-sectional planes if possible. For example, the warp-tow width, WaW, can be found by measuring warp tows from images normal to the warp direction or normal to the thickness direction. The mean values of the identified key dimensions were used in the construction of the 3D woven RUC for the model described in Section 2.3.3. In the future, automated techniques, such as those described in refs. [31][32][33][34][35], could be utilized to speed up the geometry-characterization process. Two additional features observed in the imaging of the 3D woven composite were considered in the multiscale model: disbonding of the binder tows from the surrounding matrix and waviness in the weft tows. These are depicted in Figure 3, which shows X-ray CT images of the 3D woven composite. Figure 3a is taken from a through-thickness crosssection that passes through the weft tows, and the weft-tow waviness induced by the through-thickness binder tows is evident. Figure 3b is a detail from Figure 3a, in which the disbonding of the binder tows is shown. Figure 3c is a weft cross-section showing how these binder-tow disbonds progress through the thickness of the 3D woven composite. The disbonds appear on both sides of the binder tows in the warp direction and tend to run part way through thickness from the inner radius of the binder tow as it turns from the surface of the composite to the through-thickness direction. Note that measurements of the average fiber waviness angle were taken, but it was difficult to obtain highly accurate measurements from the X-ray CT images, as the angles were quite variable, even through the width of each weft tow. Based on these measurements, the average weft-tow waviness angle was approximated to be 1.5°, to the nearest 0.5°. Note that the warp tows were far less wavy than the weft tows. Both the weft-tow waviness and the binder-tow disbonds were included in the multiscale model of the 3D woven composite, as described in Section 2.3.3, below. Two additional features observed in the imaging of the 3D woven composite were considered in the multiscale model: disbonding of the binder tows from the surrounding matrix and waviness in the weft tows. These are depicted in Figure 3, which shows Xray CT images of the 3D woven composite. Figure 3a is taken from a through-thickness cross-section that passes through the weft tows, and the weft-tow waviness induced by the through-thickness binder tows is evident. Figure 3b is a detail from Figure 3a, in which the disbonding of the binder tows is shown. Figure 3c is a weft cross-section showing how these binder-tow disbonds progress through the thickness of the 3D woven composite. The disbonds appear on both sides of the binder tows in the warp direction and tend to run part way through thickness from the inner radius of the binder tow as it turns from the surface of the composite to the through-thickness direction. Note that measurements of the average fiber waviness angle were taken, but it was difficult to obtain highly accurate measurements from the X-ray CT images, as the angles were quite variable, even through the width of each weft tow. Based on these measurements, the average weft-tow waviness angle was approximated to be 1.5 • , to the nearest 0.5 • . Note that the warp tows were far less wavy than the weft tows. Both the weft-tow waviness and the binder-tow disbonds were included in the multiscale model of the 3D woven composite, as described in Section 2.3.3, below.
The fiber volume fraction within the 3D woven-composite tows was characterized via analysis of multiple SEM images from both the warp and weft tows. The images were segmented to separate individual fibers and matrixes using standard MATLAB-based image-processing techniques (gray-scaling and thresholding).  The fiber volume fraction within the 3D woven-composite tows was characterized via analysis of multiple SEM images from both the warp and weft tows. The images were segmented to separate individual fibers and matrixes using standard MATLAB-based image-processing techniques (gray-scaling and thresholding).

Acid-Digestion Tests
Acid-digestion testing based on ASTM D3171 Procedure B [36] was used to characterize the overall fiber volume fraction and void volume fraction. A total of nine specimens were tested. The MsRM approach [6,7,25] can be used to simulate materials with multiscale microstructures. Recursive procedures, subroutines, and data structures are extensively used in the NASMAT software. This permits an arbitrary number of scales to be defined, and data can be seamlessly transferred among the scales. Additionally, any micromechanical method can be utilized at any scale. Through a series of localization operations, successively lower levels are called until, ultimately, the individual constituents are reached. Once a subvolume at any level contains a constituent material, no further localizations can occur. After all localizations have occurred at a given level, properties are homogenized and passed up to the previous level. The highest level in the model, Level 0, represents the entry point for a NASMAT analysis and is where most data of interest are generated, such as the effective composite moduli and effective stress-strain response. Currently, Level 0 calculations are performed using any of the available micromechanical theories, and periodicity is assumed. To relax the assumption of global periodicity, NASMAT can be called from a third-party tool (e.g., a finite element program) where the desired boundary conditions can be applied. Each subvolume within the Level 0 RUC can call an independent Level 1 unit cell utilizing a given micromechanical theory. This process can continue up to an arbitrary k number of levels and is schematically shown in Figure 4. Acid-digestion testing based on ASTM D3171 Procedure B [36] was used to characterize the overall fiber volume fraction and void volume fraction. A total of nine specimens were tested.

Warp and Weft Tensile and In-Plane Shear Tests
In-plane tensile and shear tests were conducted on the woven composite at the National Institute for Aviation Research (NIAR) [37]. Uniaxial tensile tests were conducted with loading in the warp and weft directions based on ASTM Standard D3039 [38], while V-notched rail shear tests were performed according to ASTM Standard D7078 [39]. Five tensile tests (per loading direction) were performed using untabbed, 254 mm × 25.4 mm rectangular specimens. Four in-plane shear tests were performed using 55.9 mm (warp direction) × 76.2 mm (weft direction) rectangular specimens with two 90°, 12.7 mm deep notches aligned in the warp direction to give a 30.5 mm width between the notch tips. Strains were measured using a strain gauge and digital image correlation (DIC).

Multiscale Recursive Micromechanics
The MsRM approach [6,7,25] can be used to simulate materials with multiscale microstructures. Recursive procedures, subroutines, and data structures are extensively used in the NASMAT software. This permits an arbitrary number of scales to be defined, and data can be seamlessly transferred among the scales. Additionally, any micromechanical method can be utilized at any scale. Through a series of localization operations, successively lower levels are called until, ultimately, the individual constituents are reached. Once a subvolume at any level contains a constituent material, no further localizations can occur. After all localizations have occurred at a given level, properties are homogenized and passed up to the previous level. The highest level in the model, Level 0, represents the entry point for a NASMAT analysis and is where most data of interest are generated, such as the effective composite moduli and effective stress-strain response. Currently, Level 0 calculations are performed using any of the available micromechanical theories, and periodicity is assumed. To relax the assumption of global periodicity, NASMAT can be called from a third-party tool (e.g., a finite element program) where the desired boundary conditions can be applied. Each subvolume within the Level 0 RUC can call an independent Level 1 unit cell utilizing a given micromechanical theory. This process can continue up to an arbitrary k number of levels and is schematically shown in Figure 4.  out of the N α i total number of subvolumes. The relation between the local and average global strains is given by: The local material's constitutive equation (assuming no inelastic or thermal effects) can be defined by: are the subvolume's stress and stiffness tensors, respectively. σ can be expressed in terms of ε i by substituting Equation (1) into Equation (2): The average (global) stress tensor, σ i , is given by: where v α i is the volume fraction for subvolume α i . Substituting Equation (3) into Equation (4) allows σ i to be expressed in terms of ε i : The effective elastic constitutive equation at Level i is given by: where C * i is the effective stiffness tensor at Level i. By comparing Equations (5) and (6), C * i is given by: In MsRM, the scales are linked by equilibrating the homogenized average stress, strain, and stiffness tensors at Level i to the local stress, strain, and stiffness tensors of a given subvolume at Level i−1 (with appropriate transformation to account for the potential coordinate system change from scale to scale). That is: where T i 2 and T i 4 are the appropriate second-and fourth-order coordinate transformation tensors, respectively. Hence, it is clear that starting with the lowest-scale (k) microstructure (see Figure 4), whose subvolumes contain only monolithic materials, the effective stiffness tensor can be calculated using any standard micromechanical theory. This stiffness tensor (after appropriate coordinate transformation) then represents the homogenized material occupying one of the subvolumes within a composite material at the next higher length scale. Given the transformed effective stiffness tensors of all subvolumes at this next higher length scale, the effective stiffness tensor of the composite at this level can be determined. This stiffness tensor can then be transformed and passed along to the next higher length scale and the process repeated until the highest length scale considered (0) is reached. As an example, for an MsRM analysis considering three length scales (0, 1, and 2), the overall effective stiffness tensor can be written using Equations (7) and (8) as: Superscripts on the bracketed terms are used to indicate that all variables within the brackets are functions of the subvolume indices from the next higher length scale. This includes lower scale volume fractions and subvolume indices. This notation is adopted to fully define the subvolume at a given scale as a function of its lower length-scale contributions. For example, the Level 2 effective stiffness tensor, from Equation (8), can be written as: where there are distinct C * 2 values for every Level 1 subvolume and distinct Level 1 composites present within each Level 0 subvolume.
In addition to multiscale homogenization, multiscale localization of the stress and strain tensors can be performed with MsRM. Multiscale localization is required to obtain local fields for handling damage (and inelasticity). For the previously described example with three length scales, the local strain tensor in an arbitrary Level 2 subvolume is expressed using Equations (1) and (8) as: Again, superscripts on the bracketed terms are used to indicate that all variables within the brackets are a function of the subvolume indices from the next higher length scale. This procedure can be repeated to determine the strain tensor for any subvolume at any length scale. The stress tensor can be found by simply using the strain tensor, along with the constitutive relationship given by Equation (2), at the appropriate length scale. Note that the MsRM implementation in NASMAT accounts for the influence of thermal and inelastic strains, but these additional effects are omitted from this section for simplicity. Since MsRM can seamlessly incorporate multiple length scales into a single analysis, it is ideal for the multiscale modeling of materials such as 3D woven composites that exhibit identifiable microstructures across multiple length scales.

Crack-Band Progressive Damage Model
While fiber filament failure was modeled as sudden (brittle) at the lowest level considered (maximum stress criteria, with failure properties given in Section 2.3.3), the crack-band model allows the matrix to be progressively damaged [27]. Herein, the crack-band model was applied at the lowest, constituent level for the resin matrix material. This version of the theory assumes that the crack band is aligned normal to the local material coordinates as opposed to being aligned normal to the direction of maximum principal stress, as in ref. [28].
The mixed-mode traction-separation law for the crack-band model follows the formulation in ref. [29]. For complete details on the current implementation, the reader is referred to ref. [7]. Quadratic failure criteria are utilized to initiate a crack band within the subvolume for a monolithic matrix material (α m ): where damage initiation is related to the normal (X m ) and shear (Y m ) matrix cohesive strengths and are the Macaulay brackets. These normal and shear strengths are numerical parameters and may not necessarily be related to those obtained from experiments. The mixed-mode traction, t M , is related to an equivalent mixed-mode strain, ε M , through a triangular mixed-mode traction-strain law, as depicted in Figure 5. ε 0 M and t 0 M are the mixed-mode strain and traction, respectively, when damage initiates, and ε f M is the mixed-mode failure strain. The orientation of the crack band can be used to determine the relationship between the tractions and the material stresses σ (α m ) [7]. The damaged , is expressed in terms of the scalar damage variable, D M . For instance, assuming that the crack band initiated normal to the x 1 direction (as governed by the first inequality in Equation (12)), the affected components of the compliance matrix would be: where E and G are the undamaged Young's and shear moduli, respectively, and ν is the Poisson ratio for an isotropic material. If damage is initiated according to one of the other inequalities in Equation (12), then the appropriate compliance matrix, strain, and stress components must be used (Equation (13)).
Polymers 2022, 14, x FOR PEER REVIEW 9 of 28 experiments. The mixed-mode traction, tM, is related to an equivalent mixed-mode strain, εM, through a triangular mixed-mode traction-strain law, as depicted in Figure 5. and are the mixed-mode strain and traction, respectively, when damage initiates, and is the mixed-mode failure strain. The orientation of the crack band can be used to determine the relationship between the tractions and the material stresses ( ) [7]. The damaged compliance matrix ( ) , where ( ) = ( ) ( ) , is expressed in terms of the scalar damage variable, DM. For instance, assuming that the crack band initiated normal to the x1 direction (as governed by the first inequality in Equation (12)), the affected components of the compliance matrix would be: where E and G are the undamaged Young's and shear moduli, respectively, and ν is the Poisson ratio for an isotropic material. If damage is initiated according to one of the other inequalities in Equation (12), then the appropriate compliance matrix, strain, and stress components must be used (Equation (13)).
Assuming the triangular mixed-mode traction-strain law, as shown in Figure 5, DM can be calculated by: The final strain energy release rate, upon complete failure, is governed by the following mixed-mode law: where GI and GII are the mode-I and mode-II-strain energy-release rates of the material, respectively. GIC and GIIC, which are input properties, are the mode-I and mode-II fracture toughnesses of the material. Since the behavior of microcracks under mode-III conditions is still an active area of research, it is assumed in Equation (15) that the matrix has no fracture resistance under mode-III cracking similar to ref. [29]. GI and GII are calculated as: Assuming the triangular mixed-mode traction-strain law, as shown in Figure 5, D M can be calculated by: The final strain energy release rate, upon complete failure, is governed by the following mixed-mode law: where G I and G II are the mode-I and mode-II-strain energy-release rates of the material, respectively. G IC and G IIC , which are input properties, are the mode-I and mode-II fracture toughnesses of the material. Since the behavior of microcracks under mode-III conditions is still an active area of research, it is assumed in Equation (15) that the matrix has no fracture resistance under mode-III cracking similar to ref. [29]. G I and G II are calculated as: where t I and t I I are the mode-I and mode-II crack-band tractions, respectively, and ε I and γ II are the mode-I (normal) and mode-II (shear) crack-band strains, respectively. G I and G II are related to the mode-specific traction-strain history and a characteristic length, l c . The total area under the mixed-mode traction-strain curve ( Figure 5) is governed by Equation (15). If the physics of the micromechanical model being employed cause the damage to localize to the size of the geometric discretization (as in finite-element analysis and HFGMC [26]), the actual discretization geometry should be used as l c . That is, the characteristic length should be set equal to the element or subcell dimension normal to the crack band. This will regulate the total energy dissipated and minimize pathological mesh-dependence [28]. The crack-band implementation within NASMAT contains an option to automatically link the characteristic length to the subcell dimensions in a consistent manner across the length scales. However, for micromechanical models that do not exhibit mesh-dependence and damage localization (such as the Mori-Tanaka method and GMC), the characteristic length may be treated as a material parameter. This latter approach has been taken herein, where the GMC micromechanical model is used at the lowest level, where the crack-band model is active. The effect of the l c parameter on the predicted 3D woven-composite behavior will be examined.

Multiscale 3D Woven-Composite Model
The multiscale 3D woven-composite RUC is depicted in Figure 6, both with and without the pure matrix regions included. While this is a fairly coarse representation of the RUC, the key features of the woven-composite microstructural geometry have been captured. These include the ratio and sizes of the warp and weft tows and the difference between the binder-tow cross-section at the surface and as it traverses through the thickness of the composite. For reference purposes, the through-thickness (TT), weft, and warp material directions are aligned with the NASMAT global x 1 , x 2 , and x 3 coordinate axes, respectively. Figure 7 shows schematically the full MsRM multiscale model for the 3D woven composite. Level 4 refers to the lowest length scale considered in the composite-that of the individual fibers and the matrix within the tows. These are homogenized to represent the tows using a 2 × 2 GMC RUC with a fiber volume fraction of 0.672 (as shown with a blue fiber and red matrix at Level 3 in Figure 7). At Level 3, the tows and the pure matrix regions between the tows are homogenized as through-thickness stacks, which is part of the double-homogenization procedure typically employed to compensate for the lack of shear-normal coupling in GMC when modeling woven composites [7]. Additionally, the appropriate fiber tow orientations are included in Level 3 RUCs. This results in a doubly periodic model at Level 1, where each through-thickness stack from Level 2 is depicted as a unique, anisotropic subcell. At this level, the HFGMC micromechanical model has been employed, which has been shown to give much more realistic in-plane shear predictions compared to GMC [25]. Finally, at the global scale, Level 0 represents the effective nonlinear response of the composite as the homogenized RUC behavior. Note that localization down the scales must also occur for every iteration at each increment of the applied loading, as damage occurs at the lowest levels (3 and 4) based on the local stresses, and the effect of the damage on the elastic properties is then homogenized up the scales.
been taken herein, where the GMC micromechanical model is used at the lowest level, where the crack-band model is active. The effect of the lc parameter on the predicted 3D woven-composite behavior will be examined.

Multiscale 3D Woven-Composite Model
The multiscale 3D woven-composite RUC is depicted in Figure 6, both with and without the pure matrix regions included. While this is a fairly coarse representation of the RUC, the key features of the woven-composite microstructural geometry have been captured. These include the ratio and sizes of the warp and weft tows and the difference between the binder-tow cross-section at the surface and as it traverses through the thickness of the composite. For reference purposes, the through-thickness (TT), weft, and warp material directions are aligned with the NASMAT global x1, x2, and x3 coordinate axes, respectively.   Figure 7). At Level 3, the tows and the pure matrix regions between the tows are homogenized as through-thickness stacks, which is part of the double-homogenization procedure typically employed to compensate for the lack of shear-normal coupling in GMC when modeling woven composites [7]. Additionally, the appropriate fiber tow orientations are included in Level 3 RUCs. This results in a doubly periodic model at Level 1, where each through-thickness stack from Level 2 is depicted as a unique, anisotropic subcell. At this level, the HFGMC micromechanical model has been employed, which has been shown to give much more realistic in-plane shear predictions The binder-tow disbonding discussed in Section 2.2.1 was incorporated along the binder tows, as shown in Figure 8. A number of thin subcells were added to either side of the through-thickness binder tows on the faces normal to the warp direction. This approximated the observed binder-tow disbonds shown in Figure 3. These disbond subcells are assigned a very low stiffness (matrix stiffness divided by 10 6 ) for all components. Note that, because these subcells have very thin dimensions in the weft direction, they would be expected to have a minimal impact on the weft-direction response.  The binder-tow disbonding discussed in Section 2.2.1 was incorporated along the binder tows, as shown in Figure 8. A number of thin subcells were added to either side of the through-thickness binder tows on the faces normal to the warp direction. This approximated the observed binder-tow disbonds shown in Figure 3. These disbond subcells are assigned a very low stiffness (matrix stiffness divided by 10 6 ) for all components. Note that, because these subcells have very thin dimensions in the weft direction, they would be expected to have a minimal impact on the weft-direction response.
As an approximation of the weft-fiber waviness discussed in Section 2.2.1, all of the weft tows were assigned an in-plane angle, with alternating columns of weft tows assigned 1.5° and −1.5°. Recall that this is roughly the average weft-fiber angle measured from the X-ray CT data. This is clearly a simplification of the actual complex weft-fiber waviness shown in Figure 3. However, it should still capture the first-order effects of the waviness, as every weft tow in the model contains an imperfect angle, which is clearly reflected in the microstructural observations. Finally, the constituent material properties employed for the RTM6 resin (treated as isotropic) and the IM7 fiber (treated as transversely isotropic) are given in Table 1. In the table, the transversely isotropic fiber properties are defined by the axial Young's modulus, E11, the transverse Young's modulus, E22, the axial shear modulus, G12, the axial Poisson ratio, ν12, and the transverse Poisson ratio, ν23. X, Y, and Z represent unique tensile, compressive, and shear strength components, respectively. Subscripts f and m are used to denote fiber and matrix properties, respectively.  Figure 6). Note that the thicknesses of these disbond subcells have been exaggerated in this figure to make the disbonds visible.   Figure 6). Note that the thicknesses of these disbond subcells have been exaggerated in this figure to make the disbonds visible.
As an approximation of the weft-fiber waviness discussed in Section 2.2.1, all of the weft tows were assigned an in-plane angle, with alternating columns of weft tows assigned 1.5 • and −1.5 • . Recall that this is roughly the average weft-fiber angle measured from the X-ray CT data. This is clearly a simplification of the actual complex weft-fiber waviness shown in Figure 3. However, it should still capture the first-order effects of the waviness, as every weft tow in the model contains an imperfect angle, which is clearly reflected in the microstructural observations. Finally, the constituent material properties employed for the RTM6 resin (treated as isotropic) and the IM7 fiber (treated as transversely isotropic) are given in Table 1. In the table, the transversely isotropic fiber properties are defined by the axial Young's modulus, E 11 , the transverse Young's modulus, E 22 , the axial shear modulus, G 12 , the axial Poisson ratio, ν 12 , and the transverse Poisson ratio, ν 23 . X, Y, and Z represent unique tensile, compressive, and shear strength components, respectively. Subscripts f and m are used to denote fiber and matrix properties, respectively.

Geometic Characterization
From acid-digestion testing, the mean overall 3D woven-composite fiber volume fraction was calculated to be 0.48 (0.01 standard deviation), and the mean void volume fraction was calculated to be 0.017 (0.003 standard deviation). Note that simulations were conducted in which diffuse matrix voids were included (although not shown herein), and it was found that such a low volume fraction of voids had a negligible effect on the composite response (see ref. [6] for details about void modeling). As such, voids were not included in the models presented in this work.
As mentioned in Section 2.2.1, a number of key dimensions were identified and measured from the X-ray CT data for use in the model RUC. These data are shown in Figure 9, with the means and standard deviations given in Table 2. Note that, aside from binder-tow measurements at the TT transition region, the measured tow widths and heights are similar. As mentioned in Section 2.2.1, a number of key dimensions were identified and measured from the X-ray CT data for use in the model RUC. These data are shown in Figure 9, with the means and standard deviations given in Table 2. Note that, aside from bindertow measurements at the TT transition region, the measured tow widths and heights are similar. Figure 9. Histograms based on measurements extracted from X-ray CT images.    Figure 10b shows a segmented image. The fiber volume fraction results are given in Table 3, and a histogram of the data is shown in Figure 11. As shown, the mean fiber volume fraction values in the warp and weft tows were quite similar, so these data were combined to arrive at a mean tow-fiber volume fraction of 0.672, which was used for all tows (including the binder) in the model of the 3D woven-composite RUC.   The dimensions associated with the geometry given in Table 2 were maintained and resulted in an overall fiber volume fraction for the 3D woven composite of 0.473 (vs. a mean value of 0.48 measured with acid digestion) when using a measured fiber volume fraction of 0.672 (see Table 3) within the tows.  The dimensions associated with the geometry given in Table 2 were maintained and resulted in an overall fiber volume fraction for the 3D woven composite of 0.473 (vs. a mean value of 0.48 measured with acid digestion) when using a measured fiber volume fraction of 0.672 (see Table 3) within the tows.  The dimensions associated with the geometry given in Table 2 were maintained and resulted in an overall fiber volume fraction for the 3D woven composite of 0.473 (vs. a mean value of 0.48 measured with acid digestion) when using a measured fiber volume fraction of 0.672 (see Table 3) within the tows.

Warp and Weft Tensile and In-Plane Shear Tests
The results from the in-plane tensile and shear tests are shown in Figure 12 and Tables 4-6. Figure 12a shows four warp and three weft tensile stress-strain curves, while Figure 12b shows three in-plane shear stress-strain curves. Some DIC data were not available near the end of the tensile tests and, as such, the plotted stress-strain curves terminated at a slightly lower stress than the measured ultimate tensile strengths (UTSs). while Figure 12b shows three in-plane shear stress-strain curves. Some DIC data were not available near the end of the tensile tests and, as such, the plotted stress-strain curves terminated at a slightly lower stress than the measured ultimate tensile strengths (UTSs). As expected, the weft-direction tensile modulus was greater than that for the warp direction because there are more weft tows than warp tows in the composite (and thus there are more fibers aligned with the weft direction). However, while the weft tensile modulus is approximately 20% higher than the warp (Tables 4 and 5), the UTSs are comparable (weft average UTS = 922 MPa, warp average UTS = 902 MPa), with a difference of only 2%.
The in-plane shear stress-strain curves shown in Figure 9b indicate a much more ductile response compared to the warp and weft tensile curves. The tests progressed to very large shear strain values, although the physical strain gauges failed earlier, while the DIC strain measurements were able to continue. These shear tests also did not result in complete failure of the specimens. Furthermore, the V-notched rail shear test standard [39] indicates that these tests are only valid up to an engineering shear strain of 5%. As such, past this point, it is likely that the results shown are not reliable as material data. As indicated in Table 6, rather than shear strength, the shear stress at an engineering shear strain level of 5% is reported. As expected, the weft-direction tensile modulus was greater than that for the warp direction because there are more weft tows than warp tows in the composite (and thus there are more fibers aligned with the weft direction). However, while the weft tensile modulus is approximately 20% higher than the warp (Tables 4 and 5 The in-plane shear stress-strain curves shown in Figure 9b indicate a much more ductile response compared to the warp and weft tensile curves. The tests progressed to very large shear strain values, although the physical strain gauges failed earlier, while the DIC strain measurements were able to continue. These shear tests also did not result in complete failure of the specimens. Furthermore, the V-notched rail shear test standard [39] indicates that these tests are only valid up to an engineering shear strain of 5%. As such, past this point, it is likely that the results shown are not reliable as material data. As indicated in Table 6, rather than shear strength, the shear stress at an engineering shear strain level of 5% is reported.

Multiscale Modeling Results and Discussion
Predictions of the nonlinear tensile and in-plane shear stress-strain responses for the 3D woven composite were made while varying the characteristic length parameter, l c . The impact of including manufacturing-induced disbonding binder tows was examined, followed by the presentation of simulations that included the observed waviness of the weft tows. Finally, the predicted responses were compared with the experimental data.
To highlight the effects of the characteristic length parameter on the response of the materials that make up the 3D woven composite, Figure 13 shows how this parameter affects the tensile and shear responses of the RTM6 matrix material, as well as the 3D woven-composite tows, which are effectively analogous to a unidirectional composite. As discussed above, a simple 2 × 2 subcell GMC RUC, with a fiber volume fraction of 0.672, was used to model the tows. Figure 13a, which displays the tensile response of the neat resin material and the transverse tensile response of the tow, shows that increasing l c leads to a more brittle response as the softening slope of the stress-strain response increases. The tow is modeled by an RUC that includes a brittle fiber (which, recall, is not subject to the crack-band model), and its response is thus much more brittle than that of the neat resin. The changing slopes in the tow stress-strain curves are due to different regions (subcells) of the matrix being damaged at different times.
In Figure 13b, in addition to the neat resin shear response (which is isotropic), both axial (σ 12 , with x 1 as the fiber direction) and transverse (σ 23 ) shear stress-strain curves are shown. As there are warp, weft, and binder tows in the 3D woven composite, the various tows will be subjected to predominantly transverse or axial shear stresses. Figure 13b shows that, again, increasing l c leads to a more brittle response, but because the applied shear loading results in mode-II-dominated simulated damage, the responses are considerably less brittle compared to the tensile loading in Figure 13a. This is because the mode-II fracture toughness is approximately four times greater than that in mode-I (see Table 1). Despite the fact that the characteristic length was varied, the mode-I and mode-II fracture toughnesses (areas under the traction-separation curves) were preserved in the crack-band model, according to the mixed-mode law in Equation (15).
It should also be noted that the responses shown in Figure 13 were generated by applying a monotonically increasing uniaxial strain. In contrast, when the resin and tow are in situ within the 3D woven-composite RUC, the loading is applied at the global scale and the local strains will be multiaxial and nonproportional during a simulation. That is, for instance, at a given point in the 3D woven composite, at a given applied global load, damage may progress significantly (e.g., local strain increasing as local stress decreases) as the model iterates. Thus, although the stress-strain curves in Figure 13 are reasonably continuous, one should not expect the global response of the 3D woven composite to be similarly free of sharp discontinuities.
in situ within the 3D woven-composite RUC, the loading is applied at the global scale and the local strains will be multiaxial and nonproportional during a simulation. That is, for instance, at a given point in the 3D woven composite, at a given applied global load, damage may progress significantly (e.g., local strain increasing as local stress decreases) as the model iterates. Thus, although the stress-strain curves in Figure 13 are reasonably continuous, one should not expect the global response of the 3D woven composite to be similarly free of sharp discontinuities.  Table 6 and different values of the characteristic length, lc. (a) Transverse tensile response. (b) Shear response.

Effect of Binder-Tow Disbonds
As described previously, pre-existing binder-tow disbonds, as shown in Figure 6, were incorporated within the multiscale model of the 3D woven composite, as shown in Figure 12, as thin subcells with very low stiffnesses. Figure 14 shows the predicted warp and weft tensile stress-strain curves, along with the predicted in-plane shear responses, for different values of the characteristic length parameter, lc, both with and without the binder-tow disbonding. From Figure 14a, it is clear that both the disbonding and the choice of lc have very minimal effects on the warp and weft tensile bahavior. This can also be seen in Tables 7 and 8, which show the predicted orthotropic effective elastic properties of the woven composite and the predicted warp-and weft-direction ultimate tensile strengths (UTSs). Eij, Gij, and νij denote individual Young's moduli, shear moduli, and Poisson ratios, respectively. Note that the effective elastic properties are independent of lc, which affects the predictions only after damage initiation (see Figure 13). As a result of the continuous tows oriented in both the warp and weft directions, the global tensile response in these directions is fairly brittle (with final failure occurring due to fiber failure). Damage does occur well before the UTS; however, its effect on the global stiffness of the composite is relatively minor.
It may be somewhat surprising that the warp direction shows such insensitivity to the binder-tow disbonding, as the disbonding is normal to the warp loading direction. However, the TT portions of the binder tows (which are disbonded) are oriented normal to the warp direction and thus do not contribute much to the composite stiffness. Furthermore, by examining the pristine warp tensile response curves in Figure 14a, it can be seen that the disbonded simulations are elastically slightly more compliant. At a stress of ap-  Table 6 and different values of the characteristic length, l c . (a) Transverse tensile response. (b) Shear response.

Effect of Binder-Tow Disbonds
As described previously, pre-existing binder-tow disbonds, as shown in Figure 6, were incorporated within the multiscale model of the 3D woven composite, as shown in Figure 12, as thin subcells with very low stiffnesses. Figure 14 shows the predicted warp and weft tensile stress-strain curves, along with the predicted in-plane shear responses, for different values of the characteristic length parameter, l c , both with and without the binder-tow disbonding. From Figure 14a, it is clear that both the disbonding and the choice of l c have very minimal effects on the warp and weft tensile bahavior. This can also be seen in Tables 7 and 8, which show the predicted orthotropic effective elastic properties of the woven composite and the predicted warp-and weft-direction ultimate tensile strengths (UTSs). E ij , G ij , and ν ij denote individual Young's moduli, shear moduli, and Poisson ratios, respectively. Note that the effective elastic properties are independent of l c , which affects the predictions only after damage initiation (see Figure 13). As a result of the continuous tows oriented in both the warp and weft directions, the global tensile response in these directions is fairly brittle (with final failure occurring due to fiber failure). Damage does occur well before the UTS; however, its effect on the global stiffness of the composite is relatively minor. tow disbond. Note that, as was the case in the experiments (see Figure 12b), the shear predictions do not exhibit complete failure or a large stress drop. Rather, the strains become quite large as the damage continues to accumulate. Note that the shear stress at 5% strain is dependent on the damage-event sequence (i.e., the jaggedness of the stress-strain curves), and, as such, the trend in the data in Table 8 is not consistent.    It may be somewhat surprising that the warp direction shows such insensitivity to the binder-tow disbonding, as the disbonding is normal to the warp loading direction. However, the TT portions of the binder tows (which are disbonded) are oriented normal to the warp direction and thus do not contribute much to the composite stiffness. Furthermore, by examining the pristine warp tensile response curves in Figure 14a, it can be seen that the disbonded simulations are elastically slightly more compliant. At a stress of approximately 400 MPa (strain of 0.007), there is a slight jog in the pristine curves that slightly reduces their stiffness and causes these curves to align closely with the disbonded binder responses. This jog is caused by damage initiated in portions of the binder tows and it appears to reduce the stiffness by a similar amount to the disbonding. In contrast to the fairly brittle predicted tensile response, as shown in Figure 14b, the 3D woven composite's predicted in-plane shear response is very ductile. Both the choice of l c and the presence of the binder disbond have a major impact on the predicted response. This can also be seen in Tables 7 and 8. The in-plane shear modulus, G 23 , is reduced by 8.1% by the binder-tow disbond, and the shear stress at 5% strain is significantly different as well. The binder disbond results in significanly lower stress after damage has been initiated, as do higher l c values, although the effect of l c is less than that of the binder-tow disbond. Note that, as was the case in the experiments (see Figure 12b), the shear predictions do not exhibit complete failure or a large stress drop. Rather, the strains become quite large as the damage continues to accumulate. Note that the shear stress at 5% strain is dependent on the damage-event sequence (i.e., the jaggedness of the stress-strain curves), and, as such, the trend in the data in Table 8 is not consistent. Figure 15 shows predicted stress-strain curves for warp, weft, and in-plane shear loading for the pristine model and models accounting for misaligned weft tows (with and without binder-tow disbonds) with increasing l c . Predicted effective properties and strengths are shown for these models in Tables 9 and 10, respectively. For warp-direction loading (Figure 15a), there is little difference between the curves, similar to the results shown in Figure 14a, and only a slight difference in the UTSs (Table 10). However, for the weft-direction loading (Figure 15b), the strength significantly decreased with increasing l c for models containing misaligned weft tows. For instance, the weft-direction UTS for misaligned weft tows without disbonds decreased from 1111 MPa to 912 MPa as l c increased from 0.008 to 0.05 (Table 10). Similar values were obtained for the misaligned weft-tow case with disbonds included as well, as shown in Table 10. In these cases, the slight misalignment appeared to make the weft tensile behavior more sensitive to the nonlinear matrix behavior, which was affected significantly by l c . The trends in the weft stress-strain curves thus followed the trends shown in Figure 13    The in-plane shear-response predictions in Figure 15 show that weft-tow misalignment has a relatively minor effect. The in-plane shear stiffness increases slightly (by 1.3%; see Table 9), and the stress-strain curves are slightly higher than the pristine curves just after damage initiation (strains between 0.015 and 0.025). In Table 10, once again, because the shear stress at 5% strain is dependent on the damage-event sequence, the trend in these data is not consistent.

Effect of Weft-Tow Waviness
When the slight weft-tow misalignment was added, overall, the effective properties were not significantly different from those in the pristine case (Table 9). However, when disbonds were added, most of the effective properties were reduced, with noticeable decreases in G 23 , G 13 , and ν 23 . This response was consistent with that of the pristine model with and without disbonds (Table 7). Figure 16 compares the in-plane tensile and shear test data with the model results for the misaligned weft tows with no disbonding, including the different values of l c . The elastic properties and strengths are compared in Tables 11 and 12. The model did an excellent job of predicting the warp and weft tensile moduli, which were within 5% of the test-data averages. The warp UTS was underpredicted by approximately 13% across all values of l c . Note that the warp UTS was observed to be insensitive not only to l c but also to the presence of binder-tow disbonds and weft-tow misalignment. As such, further investigation of both the test and model results is needed to fully understand this discrepancy. In contrast, the weft UTS sensitivity to l c in the presence of weft-tow misalignment enabled the weft UTS predictions to improve considerably with increasing l c . With l c values of 0.02 mm and 0.05 mm, the weft UTS predictions bounded the test data, with both predictions within 7%.

Conclusions
A multiscale progressive failure model for an IM7/6k carbon fiber-RTM6 epoxy 3D orthogonal woven composite has been presented. Based on the MsRM approach, the methodology employed considered the fiber and matrix constituents within the tows, as   Figure 16b shows that, while the model was able to capture the general, less brittle character of the in-plane shear response, it significantly underpredicted the in-plane shear modulus and shear stress at 5% strain. Regardless of the value of l c , the in-plane shear stress-strain curves were also underpredicted. It should be noted that the model's only nonlinear deformation mechanism is stiffness-reduction damage, whereas other nonlinear mechanisms (e.g., plasticity, viscoelasticity) could be present in shear tests. Additionally, model refinement could improve the in-plane shear predictions, particularly if the actual geometry of the weft tows (i.e., their waviness) is accurately captured.

Conclusions
A multiscale progressive failure model for an IM7/6k carbon fiber-RTM6 epoxy 3D orthogonal woven composite has been presented. Based on the MsRM approach, the methodology employed considered the fiber and matrix constituents within the tows, as well as the intertow matrix. The GMC and HFGMC micromechanical theories were used at different scales, while a crack-band progressive damage model was used for the matrix material. The model RUC dimensions and local fiber volume fractions were determined via measurements taken from X-ray CT and SEM images. Two key features observed in these images, binder-tow disbonding and weft-tow misalignment, were also incorporated into the models.
Progressive damage simulations were conducted for applied warp and weft tensions, as well as applied in-plane shear stresses. Stiffness and strength predictions, along with full stress-strain curves, were compared for different values of the crack-band length-scale parameter, l c , for cases with and without binder-tow disbonding and with and without weft-tow misalignment. It was found that both the disbonding and the choice of l c had very minimal effects on the fairly brittle warp and weft tensile behavior. However, the simulated in-plane shear response, which was much more ductile, was quite sensitive to the choice of l c and the presence of the binder-tow disbonding. The impact of weft-fiber misalignment was small on the warp-direction tensile response and relatively minor on the in-plane shear response, but quite significant on the weft-direction tensile behavior. With the addition of the weft-tow misalignment, a strong influence of the choice of l c on the simulated weft-direction tensile response was observed, indicating that the misalignment greatly enhanced the influence of the matrix nonlinearity.
Finally, comparing the predictions with experimental data, it was shown that the model's stiffness predictions were in very good agreement with the test data for in-plane tension, while the in-plane shear modulus was underpredicted to some extent. Tensile strength predictions were also reasonable, with the weft-tow misalignment enabling muchimproved strength predictions for the weft direction. Most notable was the clear ability of the presented multiscale model to capture the fairly brittle response for in-plane tension, while also capturing the much more ductile in-plane shear response observed in the experiments.  Data Availability Statement: Some data presented in this study are available on request from the corresponding author. Some data are not publicly available due to NASA and U.S. government restrictions.

Conflicts of Interest:
The authors declare no conflict of interest.