Abstract
This article presents the development of a finite element model, which considers stochastic tow waviness using a Markov Chain algorithm and non-linear material properties using Binary Model, to predict the stress–strain and fracture behaviour of plain-weave ceramic matrix composites under uniaxial extension. The stochastic waviness is described by fluctuations in the centroid coordinates of tow positioning. The tow deviations are generated by marching sequentially from one grid point to next along a tow path. The deviations depend only on the deviation of the previous point using a probability transition matrix. A non-linear orthotropic constitutive model was implemented in a commercial finite element code Abaqus using a user-defined subroutine. Two 2 × 2 unit cell models of a plain-weave ceramic matrix composite laminate are created using stochastic tow elements generated by the virtual specimen generator, which was developed on the basis of the Markov Chain algorithm. A comparison has been made between the systematic and stochastic models to assess the effects of stochastic tow waviness on the stiffness and strength of the laminate. The numerical results have been validated by the comparison of predictions with the experimental data. The stochastic model which considers random waviness correlates well with the experimental data.
Keywords
Introduction
Ceramic matrix composites (CMCs) are characterised by low density, high modulus and good thermal stability; they are of increasing interest for hot structures, such as rocket nozzles, combustor liners of turbine engines, space shuttle thermal protection systems and nuclear fuel cladding tubes.1,2 For actual fibre tows, there are some geometrical defects like matrix voids, fibre waviness and variations in tow positioning and cross-section, which may arise during the manufacturing process. Fibre waviness is one of the most frequently encountered process-induced defects in CMCs and seriously affect the performance of CMC components under complex loading conditions. Therefore, it necessitates the development of an accurate yet highly computational efficient model, which is capable of stochastically modelling fibre waviness.
In the last three decades, numerous efforts have been made to investigate the effects of fibre waviness on mechanical behaviour of woven composites. Yurgatis. 3 and Jelf and Fleck 4 measured the local fibre misalignment angle and showed that the geometric morphology of wavelength is sinusoidal. Hsiao and Daniel 5 proposed an analytical model to investigate the compressive behaviour of thick composites with fibre waviness. Chun et al. 6 developed analytical model, in which the geometric non-linearities triggered by out-of-plane waviness were considered and the tensile and compressive behaviour of composites was predicted. Zhang and Hayhurst7–9 developed a conceptually simple and computationally efficient finite element model to study the effects of fibre misalignment of 7° on stiffness and strength of CMCs under uniaxial tension. However, these works deal with fibre waviness by the means of averaging tow orientations, without considering random irregularity in realistic materials. Stochastic waviness can initiate local stress concentration and fatigue damage initiation which has significant influence on the life prediction of CMC components that would experience multiaxial stress states. Therefore, stochastic distribution of fibre waviness must be taken into account and explicitly modelled to predict accurate stress distributions and ultimate strength.
Generally, existing approaches for representing the stochastic geometry can be categorised into two types: non-intrusive and intrusive. Non-intrusive methods use the deterministic representations of geometry, while intrusive techniques involve the reformulation of model. 10 The most common non-intrusive model is in the Monte Carlo framework, which is used to generate random samples of input variables from statistical distributions. In the case of intrusive techniques, such as that of Sasikumar et al., 11 they used Karhunen–Loeve (K-L) and polynomial chaos expansions to represent the input and output random fields, respectively, to quantify the random variables of multilayer composites. In an analogous approach, Wang and Wang 12 derived the theoretical model by introducing Taylor’s expansion of local stochastic variable into stiffness of textile composites. Nevertheless, intrusive techniques require a quite large number of equations especially in the case of complex and multi-directional stochastic problems, leading to high computational costs.
This work utilises a probabilistic collocation method, that is, the non-intrusive Markov Chain algorithm, in which stochastic waviness can be described by fluctuations in the centroid coordinates of tow positioning, and then, the tow deviations are generated by marching sequentially from one grid point to next along a tow path. The deviations depend only on the deviation of the previous point, using the probability transition matrix (PTM) that can be calibrated using experimental data. Bale et al., 13 Blacklock et al. 14 and Rinaldi et al. 15 indicated that the Markov Chain algorithm is physically appropriate for reconstruction problem of textile composites. The dominant correlations must be those along a tow, while correlations between inter-tows are relatively weak, which satisfies the Markovian procedure. Consequently, the one-dimensional (1D) representations of tows with stochastic waviness are generated within virtual specimen generator.
As shown in Figure 1, the Binary Model consists of two virtual components, 1D tow elements and three-dimensional (3D) effective medium elements, which represent the fibre-dominated and the matrix-dominated properties, respectively. Almost all works using the Binary Model assumed that effective medium elements are taken to be isotropic and tow elements are considered to be linear elastic. Such as those of Flores et al., 16 Rajan et al. 17 and Rossol et al.; 18 they used a linear Drucker–Prager yield criterion for failure initiation in the effective medium along with isotropic elastic–plastic evolution law. Yang and Cox 19 employed the rule of mixture to determine the material properties and assumed that the effective medium is isotropic and the properties of both virtual constituents are linear elastic. The assumption of isotropic medium is reasonable for polymer matrix composites (PMCs) since tow elements usually take dominant role due to the low stiffness of polymer matrix. Nevertheless, this postulation may not be appropriate for CMCs because of the high stiffness of ceramic matrices and their progressive damage behaviour, which can be better described by orthotropic effective medium properties. The axial response of elastic fibres is characterised by the axial stiffness of the tow elements. The strength of fibres can be taken to be either deterministic or stochastic. Most work using the Binary Model assumed that the properties of tow elements are linear elastic using the simple choice of deterministic strength. If the fibre strength is stochastic, the cumulative effects of fibre-matrix debonding, fibre fractures and subsequent frictional pull-out within each tow are stochastically modelled by a 1D non-linear constitutive law. In this article, the latter method has been used. On the basis of previous work,7–9 this article employs a non-linear orthotropic constitutive model. Consequently, accurate predictions of their mechanical behaviour, for example, non-linear stress–strain curves, orthotropic damage initiation and evolution and catastrophic failure in a brittle manner can be obtained.

Principle of binary model: one tow element AB embedded in one eight-noded effective medium element.
This article addresses the development of a finite element model which considers both the stochastic tow waviness and the non-linear tow properties. This work relies on a good understanding of characteristic of tow misalignments as well as constitutive properties of the individual composite constituents. The layout of this article is as follows. First, the microscopic structure of CMCs will be introduced. Second, the statistical analysis of stochastic waviness using the Markov Chain algorithm is made. Then, the Binary Model and its non-linear orthotropic constitutive equations employed are described. After that, the test data of the 10 high-carbon fibre/amorphous carbon matrix–SiC matrix (C/C-SiC) DLR-XT plain-weave laminates will be used to validate the predictions. Also, the effects of stochastic waviness will be investigated on composite stiffness, strength and overall stress–strain response.
Material description
The microstructures of CMCs depart significantly from those of PMCs due to their constituents and infiltration method. Polymer matrix can be adequately infiltrated into interior tows and interstices of fibre bundles, which creates nearly porosity-free and homogeneous matrices. For CMCs, in contrast, the matrix is formed by chemical vapour infiltration (CVI) or other process that commonly contains residual microcracks and porosities over relatively large volumes, as illustrated in Figure 2. The distribution of matrix within and between the fibre tows appears to be non-uniform. Additionally, Figure 2 shows the wavy tow paths, which are idealised by piecewise linear tow elements. These paths of a single layer are marked by yellow broken lines.

Side elevation of microstructures for DLR-XT material required from SEM images, along with idealised tow topologies of a single layer.
Statistical analysis of stochastic waviness
The reconstruction algorithms or generators for stochastic heterogeneous materials are related to the statistical physics. 14 A statistical description of tow centroid deviations in the plain-weave architecture is made. The statistical parameters required to calibrate the virtual specimen generator are introduced in this section. The formulation of Markov Chain algorithm is also presented.
Systematic variations and stochastic deviations
The unit cell of plain-weave DLR-XT material consists of two warp and two weft yarns with reference periods

Schematic drawing of the unit cell with reference periods for DLR-XT material, all dimensions in millimetre.
A uniform grid is defined over the reference period with nine grid points, i = 1, …, 9. It is assumed that the tow at grid point 10 is statistically equivalent to the tow at point 1. All tows are discretised with an equal number of grid points and the grid spacing of tow genus. The coordinates of the centroid of any section of individual tow can be represented by the vector
where
Due to the periodicity of systematic variation, any warp or weft tow can be made to coincide with any other tow of the same genus by shifting some vector (
Determination of the statistical parameters
The deviations of centroid coordinates from the systematic value at all available points can be used to (1) determine the root mean square deviations (RMSDs) and (2) determine the correlation lengths within one tow and inter-tows.
The RMSD
with
Correlation length is defined in terms of Pearson’s correlation parameter for pairs of data
where n is the number of pairs. The correlation length is essential to evaluate the spatial dependencies that means deviations at a certain position have an influence on the deviations at another position. For each tow genus, this parameter is determined by linear approximation of the first autocorrelation values for small spacing of grid points
with
The RMSD and correlation length parameters for the various components of the deviations of the tow centroid coordinates are summarised in Table 1.
Experimental data: RMSD and correlation length for the in-plane and out-of-plane deviations of centroid coordinates of tows.
Generation of tow centroid deviations with the Markov Chain algorithm
The tow centroid fluctuations around the mean centroid paths are generated from the Markov Chain algorithm. As for the angle interlock weave of Blacklock et al.,
14
no significant inter-tow correlations for the same genus are presented. Hence, the plain-weave composites do not exhibit significant inter-tow correlations. The deviations of considered parameter are discretised on an interval { –ma, −(m − 1)a, …, 0, …, (m − 1)a, ma} with grid spacing a and number of intervals m that satisfy the relation ma = 3σ. The possible values of generated deviations are limited by the 3σ, since the high-standard deviation may cause the high amplitude of spikes in the generated deviations along the tow length. The probability of occurrence of the discrete values of
where T denotes the transpose operation. The Markov process generates the distribution vector
In the initial form, this PTM is constructed as (2m + 1) × (2m + 1) tridiagonal PTM
All elements have values between [0,1] and each column adds to unity. The generated deviations
The Markovian process assumes that a deviation at any point only depends on the deviation at the previous point and not on the distribution of the deviation at prior points. Since the PTM is assumed to be independent for tows in a given genus, the process is stationary. To eliminate grid effects arising from the Markov Chain operator, the stochastic deviation can be smoothed to remove mesh-dependent small-scale fluctuations using the modified moving average.
Finite element model
In order to demonstrate the fidelity of virtual specimen generator to create large numbers of replicas of the microstructures, two samples of 2 × 2 unit cell Binary models are generated for the plain-weave DLR-XT laminate. The material properties of effective medium and tow elements are the required input for Binary Model simulations. The effects of the random waviness on the stiffness and ultimate strength of the laminate are studied. Finally, a comparison is made between the systematic and stochastic predictions as well as the results of more conventional finite element models and the experimental data.
Constituent properties for the Binary Model
The Binary Model utilises two virtual components that comprise the composite tow: a 1D tow element and a 3D effective medium element to represent the fibre-dominated and the matrix-dominated properties, respectively. 20 Effective medium elements, which represent matrix-dominated contributions to stiffness, are assumed orthotropic with nine independent material properties. Due to the reinforcement of continuous fibre bundles, CMC tows exhibit significantly different longitudinal and transverse non-linear behaviour. The degradation of different longitudinal and transverse shear stiffness as well as Poisson’s ratios can be taken into account. Tow elements represent the fibre-dominated properties and therefore are assumed to have only axial stiffness. Unlike the use of deterministic strength of elastic fibres, the cumulative effects of fibre-matrix debonding, fibre fractures and consequent frictional pull-out mechanisms within each tow can be stochastically modelled by a non-linear constitutive law.
The failure criteria for effective medium and tow elements are defined by the maximum principal strain criterion to capture the strain-induced damage, which means when strain components in principle material direction exceed the critical value associated to the cracking, the damage takes place.
The constituent properties used in the Binary Model are derived from the material properties for a DLR-XT unidirectional lamina reported in Zhang and Hayhurst.
7
In the local coordinate system (detailed in Appendix 1), equation (9) implies that there are many possible solutions for the axial stiffness of a tow element
Substitution of
It can then be determined that

Material properties of effective medium for DLR-XT composites: (a) longitudinal stress–strain curve; (b) Poisson’s ratios,

Material properties of tow elements for DLR-XT composites.
The non-linear constitutive properties of tow elements and orthotropic constitutive properties of effective medium elements are implemented by Abaqus/standard 21 with a user-defined field subroutine (USDFLD), in which field variables are used to govern the degradation of constitutive properties.
2 × 2 unit cell Binary Model generated from virtual specimen generator
In comparison with other numerical models for woven CMCs, for example, conventional fine-mesh model 22 or homogenised coarse-mesh model, 7 the mesh generation of Binary Model is much more straight forward. The topology of woven tows can be represented by 1D tow elements. Due to the periodic and symmetrical characteristics of plain-weave DLR-XT laminate, a unit cell has been chosen to represent the whole laminate. The Binary Model employs two mesh systems comprising two-node truss elements representing the fibre tows and eight-node solid effective medium elements that define the external geometry and represent the matrix-dominated properties. The two mesh systems facilitate the spatial matching of fibre tows and matrix and obviate the complexity of matrix grids, which are coupled by the multi-point constraints in Abaqus. To satisfy the continuity of stress and compatibility of displacement at the boundaries of unit cell, periodic boundary conditions have been applied.
In addition, the stochastic 1D tow representations that are generated from virtual specimen generator are suitable for insertion into Binary Model simulations of textile composites.14,23 On the basis of statistical parameters defined in Table 1, two 2 × 2 unit cells are generated for the plain-weave DLR-XT laminate. The difference between the two finite element models lies in the tow waviness, where the systematic model has no waviness and the stochastic model has random waviness, as shown in Figures 6 and 7.

Mesh and geometry of a systematic 2 × 2 unit cells Binary Model of a plain-weave DLR-XT (C-C/SiC) material, all dimensions in millimetre.

Mesh and geometry of a stochastic 2 × 2 unit cells Binary Model of a plain-weave DLR-XT (C-C/SiC) material, all dimensions in millimetre.
In the solution of a non-linear finite element analysis, the non-linear stress–strain properties were discretised to multi-linear curves. The loading was imposed in terms of the displacement boundary condition, that is, loading was modelled with displacement control rather than force control. The applied displacement was divided into many small increments. In each increment, linear material properties were used in the constitutive equations for both tow element material and for effective medium/matrix material (detailed in Appendix 2).
Results and discussion
Fibre waviness exists in woven CMCs, and in fact, it is often significantly more severe than that in unidirectional composites, especially for high fibre volume fraction. Figure 8 shows the stochastic in-plane and out-of-plane centroid coordinate paths around the systematic paths for warp genus. The out-of-plane centroid coordinates appear to have the largest deviations, as observed in Figure 8(b). This is due to the compression and interaction of preform through its thickness during the weaving process, and a tow can rotate more easily lying out-of-plane than in-plane. The relatively straight tow segments do not exhibit substantial misalignments, which might be attributed to the presence of flat surface in the moulding process.

Stochastic deviations from the systematic path for the centroid coordinates of the warp tows in 2 × 2 unit cells Binary Model of a plain-weave DLR-XT (C-C/SiC) material: (a) in-plane centroid paths and (b) out-of-plane centroid paths.
It can be seen that all of the tow centroid coordinate paths vary in the deviations, rarely exhibiting the symmetry implied by their representations in the systematic path. All deviations are less than
Using the developed Binary Model and the non-linear orthotropic constitutive properties shown in Figures 4 and 5, two sets of predictions are made for the plain-weave DLR-XT laminate under strain-controlled uniaxial tension. Figure 9 shows a comparison among the systematic and stochastic predictions, the finite element results 7 and the experimental data in the literature. 24 It can be seen that the effects of waviness on the stiffness of the stress–strain curves are significant. The stochastic model predicts a lower stiffness, a lower ultimate strength, but a higher failure strain than the systematic model. Upon loading, both curves predicted by systematic model and stochastic model are similar until the proportional limit is reached. But afterwards, the curve for the stochastic model yields significantly lower stiffness and stress. The range of tensile strengths predicted is from 210.27 to 190.26 MPa, which indicates a difference of approximately 10%. This is because the well-aligned fibre tows of systematic model can carry greater stresses than misaligned tows, and therefore, the former model predicts an earlier failure, but a higher strength than the latter one.

In comparison with experimental data, only the prediction which considers the waviness compares well with experimental data. For both the systematic and the stochastic models, the predicted failure strains have good correlations with the test data. The catastrophic manner observed in the experiments is also well captured by the sharp drop in the stress at the failure strain. The same conclusions have been drawn in Zhang and Hayhurst. 7
Conclusion
Two binary models of 2 × 2 unit cells have been created by the virtual specimen generator that possesses the same statistics as experimental data derived from scanning electron microscope (SEM) images. The Markov Chain algorithm with specific input parameters generates deviations of the coordinates of tow centroids that subsequently are added to systematic path. A stochastic model with random waviness is obtained. Comparisons between systematic and stochastic models have been made to investigate the effects of random waviness on stiffness and strength of the stiffness and ultimate strength of CMC laminate under uniaxial tension. The accuracy and computational efficiency of the virtual specimen generator have been evaluated by comparing the predicted and experimental results. Some conclusions can be obtained as follows:
The non-linear orthotropic Binary Model considering stochastic waviness can accurately predict the axial tensile behaviour of a plain-weave DLR-XT laminate.
The out-of-plane centroid coordinates are subject to the largest deviations from its systematic path due to the compression and interaction of preform through its thickness during the weaving process.
The random waviness has significant effects on the material stiffness and strength. The loss of tensile strength is approximately 10%.
The high efficiency of reconstructing the same statistics of measured sample makes the approach a useful tool to model large multi-scale-woven composites.
Footnotes
Appendix 1
Appendix 2
Academic Editor: Farzad Ebrahimi
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (11272207) and the Specialized Research Fund for the Doctoral Program of Higher Education (20120073120019).
