Abstract
Microstructural fracture behavior of a ceramic matrix composite (CMC) with nonuniformly distributed fibers is studied in the presentation. A comprehensive numerical analysis package to study the effect of nonuniform fiber dimensions and locations on the microstructural fracture behavior is developed. The package starts with an optimization algorithm for generating representative volume element (RVE) models that are statistically equivalent to experimental measurements. Experimentally measured statistical data are used as constraints while the optimization algorithm is running. Virtual springs are utilized between any adjacent fibers to nonuniformly distribute the coated fibers in the RVE model. The virtual spring with the optimization algorithm can efficiently generate multiple RVEs that are statistically identical to each other. Smeared crack approach (SCA) is implemented to consider the fracture behavior of the CMC material in a mesh-objective manner. The RVEs are subjected to tension as well as the shear loading conditions. SCA is capable of predicting different fracture patterns, uniquely defined by not only the fiber arrangement but also the specific loading type. In addition, global stress-strain curves show that the microstructural fracture behavior of the RVEs is highly dependent on the fiber distributions.
1. Introduction
Ceramic materials have been used in high-temperature applications due to their excellent thermal properties. However, the ceramics are highly brittle in nature and can store little amount of strain energy when they are stressed. In order to increase toughness of the ceramic materials, they are often reinforced by ceramic or carbon fibers, which are commonly known as ceramic matrix composites (CMCs). Since the CMC materials still exhibit superior thermal stability with enhanced toughness, they are utilized for structures exposed to high-temperature environments, e.g., blades of a gas turbine, thermal protection systems of a space vehicle, and aircraft disc brakes [1]. The improved toughness of the CMC materials is mainly due to the existence of interfacial coating layers between the fibers and matrix material. In general, the reinforcing fibers are coated with various materials to protect them from thermal oxidation and preserve chemical stability at high temperatures. The coating layer is essentially a compliant interphase. It can arrest cracks, deflect crack propagation paths, and prevent cracks occurring in the fiber and matrix phases from joining together [2]. As a result, the complicated fracture process significantly improves the damage tolerance of the CMC materials although each constituent is still brittle [3, 4]. The complex fracture behavior occurring inside the material at the fiber level is typically manifested as a nonlinear stress-strain response at the macrolength scale. Therefore, it is crucial to have a better understanding of the microscale fracture behavior to predict the global mechanical response of the CMC materials.
Micromechanical analysis of composite materials reinforced by unidirectional continuous fibers is typically performed either using analytical schemes or using representative volume element (RVE) models. Historically, Hill [5] and Hashin [6] developed analytical homogenization methods to compute elastic constants of the composites by constructing a concentric cylinder array and solving a boundary value problem. Homogenized transverse shear modulus was not singly determined, but its upper and lower bounds were calculated from their approaches. Christensen and Lo [7] later obtained an approximate solution for the transverse shear modulus using Eshelby’s inclusion theory, which tracks the strain energy change due to the existence of an inclusion in an infinite medium [8]. Suctu [9] presented a recursive algorithm to calculate homogenized elastic constants of the composites reinforced by unidirectional fibers coated with arbitrary number of layers. Guinovart-díaz et al. [10] developed a recursive asymptotic homogenization scheme (RAHS) to estimate the effective properties of the composite materials more precisely. The predictive capability of RAHS has been demonstrated in [10], but the scheme is limited to solve well-defined boundary value problems with a regular array of fibers.
RVE models defined in the transverse cross section of the composite materials have successfully been implemented into the numerical analysis for predicting the effective elastic properties of the fiber-reinforced composites. Sun and Vaidya [11] were the first to use the RVE approach to compute effective composite properties. They considered RVE models with simple square and hexagonal fiber arrays. Since then, in order to consider nonuniform or realistic fiber distributions in the RVE models, various algorithms have been developed by many researchers. Many algorithms are based on the random sequential adsorption (RSA) model or the hard-core model (HCM) [12]. Yang et al. [13] used transverse section images of a CMC material and measured statistical distributions of fiber diameters and locations. They utilized HCM in conjunction with the measured statistic data to create an RVE for the CMC material. Merlo et al. [14] also utilized HCM but relocated the initial fiber positions using their “stirring” algorithm to create resin-rich areas and achieve higher fiber volume fractions. Yang et al. [15] developed the random sequential expansion (RSE) algorithm that is based on the hard core model but overcomes the jamming limitation of HCM. Fibers are placed around the first fiber until no more space is available for placing an additional fiber. The distances between fibers can be determined based on experimental measurements or empirical estimation.
There also exist various algorithms that are based on logics other than HCM. Vaughan and McCarthy [16] utilized the nearest neighbor algorithm (NNA) to fill the RVE area with fibers sequentially with the constraint on the distances between fibers. The distance is randomly selected from the distribution of the statistical data of the microstructure measured from a section image. Wang et al. [17] further refined NNA by adding an additional criterion before removing a fiber that does not satisfy the nearest distribution. Ismail et al. [18] first created structured arrays of fibers and then changed fiber diameters and their location to create an RVE with a higher fiber volume fraction. The fiber diameters are randomly selected from the probability density function determined from experimental measurements. The fiber locations are altered by imposing velocities on each fiber and the collision between fibers are modeled with the Hertz contact law. Weng et al. [19] also started from a structured hexagonal fiber array and assigned velocities to the fibers. The elastic collision model is applied to create randomness in the fiber distribution.
In addition to the 2D RVE models for unidirectionally continuous fiber reinforced composites, the RVE modeling approach is also applied to construct 3D models having reinforcements of various configurations. The random sequential adsorption (RSA) algorithm has also been applied to generate 3D RVE models with various types of reinforcements such as short/long fibers and ellipsoidal inclusions [20–27]. Pathan et al. [28] started from the hard-core model but utilized an optimization technique to remove fiber jamming. They defined a penalty function with the total fiber intersection and minimized the function by controlling fiber coordinates until nonintersecting fiber arrangement was achieved. Catalanotti [29] developed an algorithm to obtain a high-volume fraction composite by perturbing a regular packing of particles until the statistical criteria were fulfilled. Bahmani et al. [30] utilized the event-driven molecular dynamic (EDMD) approach. By mimicking the collision of molecules, they generated high fiber volume fraction unidirectional composites with the realistic distribution of fibers.
Aforementioned numerical models are all concerned with computing effective elastic properties through the RVEs with nonuniform fiber distributions yet following the geometric information of the real microstructure. Such RVEs have been utilized to study the effects of the nonuniformity of the fiber arrangement on the stress distribution [31, 32]. Hojo et al. [33] studied the effect of the spatial arrangement of fibers on the interfacial normal stresses between fibers for thermally and transversely loaded fiber-reinforced epoxy materials. They found the maximum interfacial normal stress at the location where the distance of the fiber is the shortest and the angle between them is aligned with the loading direction. Romanowicz [34] studied microscale progressive failure behavior of a unidirectionally fiber-reinforced polymer composite subjected to transverse tensile loading. A cohesive zone model with a bilinear-type cohesive law for mixed mode debonding was utilized to model the fracture behavior. Additionally, he considered interphase layers, of which properties are radially graded. Only one case of an RVE having randomly distributed fibers was considered in his study. Merlo et al. [35] studied the microstructural fracture behavior of polymer composites reinforced by unidirectional fibers with various loading conditions. They considered the plastic deformation of the matrix material and debonding at the fiber/matrix interface. Their RVE models had random fiber distributions, but they were not connected with any experimental measurements. Meyer and Waas [2] performed micromechanical analysis of CMC materials subjected to transverse tension. They considered fracture behavior of the matrix and coating materials using the crack band method. Two baseline RVE models were created based on section images, and the fiber sizes and locations were altered to study the effect of the randomness of the fiber size and location. However, any statistical relevance of the random models to the baseline models or experimental observations was not reported.
In addition to the previous works mentioned here, there exist numerous examples of the RVE-based micromechanical analysis in the literature. Similarly, with the studies listed here, they have been mainly focused on the development of a technique that can generate realistic or statistically equivalent nonuniform fiber distribution and the computation of effective elastic properties. Beyond the continuum regime, failure initiation and consequent microscale fracture behaviors dependent on fiber arrays are not yet broadly covered. Furthermore, previous studies on the microstructural failure behavior with the nonuniform fiber arrangement are mainly concerned with polymer composites reinforced with unidirectional continuous carbon or glass fibers. There are few studies concerning the CMC materials of which the microstructure is more complicated than the traditional polymer-based composite materials due to the additional coating layer. Here, we present the study on the effects of the spatial arrangement of the fibers on the microstructural fracture behavior occurring inside the CMC materials subjected to various loading configurations. In doing so, we have developed a new algorithm to generate nonuniform fiber distributions conforming to statistical constraints determined based on the experimental measurements.
2. Characterization of a CMC Microstructure
The present numerical package starts from constructing a square RVE with the side length of H that contains the randomly distributed fibers of different dimensions as shown in Figure 1. The fibers in the RVE are modeled as a circle, and the coatings are assumed to be concentric with the corresponding fibers. Such an RVE is supposed to represent an actual microstructure of a CMC material as shown in Figure 3 of [1]. In order for the RVE to be the representative of the CMC microstructure, the volume fractions of each phase in the RVE should first comply with those of the material. In the presentation, the volume fractions of the fibers and coatings are measured from the binary image as shown in Figure 2. The binary image in Figure 2 results from the SEM image in Figure 3 of [1]. The coating layers in Figure 2 are filled with white color while the fibers and the matrix are displayed in black color. The numbers in Figure 2 are the pixel numbers of the individual fibers and coatings, which can be converted to the areas of each phase. From the areas, the volume fractions of the fibers and the coating layers in Figure 2 are found to be 0.35 and 0.06, respectively. The physical areas are also used to compute the diameters for the corresponding circular fibers and the coatings for the RVE model, which result in the same volume fractions of each phase. Figure 3 shows the computation results for the equivalent fiber diameters and the coating thicknesses obtained from the actual CMC microstructure in Figure 2. The measurement results for the fiber diameters and the coating thicknesses can be summarized as 13.527 ± 1.053 μm and 0.800 ± 0.221 μm, respectively. The mean values are also marked in Figure 3.
3. RVE Construction Process
3.1. Random Generation of the Fiber and Coating Dimensions
In this section, the procedure for constructing an RVE with random fiber dimensions and locations will be described in detail. The measurement data in Figure 3 will be utilized for demonstrating the numerical process. The RVE modeling approach is typically employed mainly because computational costs can be significantly reduced without losing accuracy in micromechanical analysis. The computational efficiency is achieved through the reduction of the model size and the accuracy is assured when the RVE model preserves key features of a microstructure. One of the important parameters is the volume fractions of constituents.
The present numerical framework also starts from reducing the model size. The reduced number of fibers to be included in the RVE model is a user-specified value. Once the fiber numbers are determined, the fiber diameters are randomly generated by a random number generator (RNG) within the range bounded by the minimum and maximum measurements. However, since the random number generation is solely based on the probabilities of each random value in the particular range, the resultant dimensions would not perfectly duplicate the statistical distribution as well as the volume fraction of the original microstructure. Therefore, the random fiber diameters need to be calibrated to represent the actual microstructure. Since the fiber dimensions should simultaneously satisfy multiple constraints including the mean, the standard deviation, and the fiber volume fraction, the calibration process can be transformed into an optimization problem.
Here, the problem is set to find an optimal random set of fiber diameters that produces the identical fiber volume fraction measured from the actual microstructure. The zero tolerance is strictly enforced on the fiber volume fraction because the volume fraction predominantly determines the representativeness of the reduced model. Furthermore, the volume fraction acquired from the image in Figure 2 is a determinate value without any deviation since no simplification on the geometries is applied in computing the areas of each constituent. For a square-shape RVE model with the side length of as illustrated in Figure 1, the fiber volume fraction, , can be computed fromwhere is the vector composed of , which is the randomly generated fiber radius. denotes the fiber index, and is the number of the fibers in the RVE model. is generated in the range between and , the minimum and the maximum measurements, respectively. When is the measured fiber volume fraction of the original microstructure, the optimization problem can be formulated aswith the objective function being defined as a sum of squared errors for the mean and the standard deviation:where and are the mean and the standard deviation of the design variables, respectively. and are, respectively, the mean and the standard deviation of the fiber measurements in Figure 3. The subscript denotes quantities associated with the fiber. The objective function, , is always positive except when the mean and the standard deviation are exactly the same as those of the real microstructure. In addition to the volume fraction constraint, the optimization formulation in equation (2) restrains the errors of the mean and the standard deviation, and , respectively, while the minimization process is being performed. and are defined as
If an optimal solution that satisfies all the constraints is not found from the first attempt, the optimization process starts over from the generation of a new set of random numbers.
The coating dimensions for the RVE model are obtained in a similar fashion. When the subscript denotes quantities associated with the coating and is the vector composed of the coating thicknesses, the optimization formulation can be expressed aswhere the objective function is defined asand the errors for the mean and the standard deviation are
The coating volume fraction in the reduced model can be computed fromwhere is the randomly generated coating thickness. , , and are the fiber radius, the fiber index, and the number of the fibers in the RVE model, respectively, as already defined in equation (1). The fiber diameters are assumed to be already determined from equation (2) before the coating thickness search is performed, and thus the fiber radius is not a design variable of the optimization process in equation (5). Since the minimization problems (equations (2) and (5)) include both equality and inequality constraints, they are solved using the interior point method [36].
3.2. Random Distribution of the Coated Fibers
After the dimensions of the coated fibers are fully determined, they are randomly placed in the RVE area. The center coordinates of each fiber are randomly generated using a uniform distribution random number generator. Since the random placement does not account for any interaction between the fibers, some fibers may intersect with others as illustrated in Figure 4(a). Such a fiber intersection, however, is physically impossible to exist in reality. In order to remove the unrealistic intersecting area without disturbing the randomness of the fiber locations, any intersecting pair is assumed to be connected through a virtual spring as shown in Figure 4(a). The virtual spring of the two coupled and fibers in Figure 4(a) is considered to be in a compressed state and thus produces a virtual repulsive force (Figure 4(b)) that eventually separates the two fibers.
(a)
(b)
Figure 4(b) shows the virtual spring with the reference coordinate system of -axis and -axis and the local coordinate system of -axis and -axis. The coordinate system is rotated by an angle of about the origin of the coordinate system. At the center points of the and fibers, the horizontal and vertical degrees of freedom (DOFs) are designated with and , respectively, according to the global coordinate system. The subscripts and of the DOFs denote the and fibers, respectively. In a similar manner, the DOFs, and , with respect to the local coordinate system are also assigned at each node.
The virtual repulsive forces at each fiber center are illustrated in Figure 4(b). The virtual force, , of the fiber with respect to the coordinate system can be expressed aswhere is the spring constant or can be considered as a dispersion coefficient. is the distance between the two fibers, is the radius of the fiber, and is the radius of the fiber as shown in Figure 4(a). is a small supplemental displacement to ensure complete separation between the two fibers. If the virtual spring is in the state of equilibrium, the force at the fiber can be easily obtained as
Forces in the transverse direction, and , are zeros under the local coordinate system. When the virtual spring is considered as a one-dimensional rod element in the typical finite element framework, the analogous local matrix system for the spring element shown in Figure 4(b) can be constructed asor in a shortened form,where is the force vector, is the local element stiffness matrix, and is the displacement vector with respect to the local coordinate system. The displacement and the force vectors in equation (11) can be related to the corresponding quantities in the global coordinate system through the angle of , and it reads that
Equations (13) and (14) can be expressed in a shortened form aswhere is the displacement vector and is the force vector under the global coordinate system. is the transformation matrix. Substituting equations (15) and (16) into equation (12) and multiplying both sides by the inverse of the transformation matrix, , yieldswhere the stiffness matrix in the global coordinate system, , is expressed as
Knowing that the transverse forces, and , are zeros under the local coordinate system, and using equation (10), the simplified expression of the force vector with respect to the global coordinate system can be obtained after multiplying both sides of equation (14) by the inverse of the transformation matrix:
Equation (17) exists for any intersecting fiber pairs. Assembling all such local matrices result in the global matrix system that can be written aswhere is the global force vector, is the global stiffness matrix, and is the global displacement vector. Equation (20) is first solved for . Then, and are updated since in equation (9) and in equation (18) have been changed due to the displacement of fibers obtained from . This process is repeated until the L2-norm of in equation (20) is sufficiently close to zero, which means all the fibers are separated without any interference.
While the fibers are not allowed to intersect with each other, intersections between coating layers may exist as observed from the section image in Figure 2. When this happens in the RVE model after the fibers are relocated due to the virtual springs, the total area of the coating layers would not be the same as that of the optimal dimension, , computed from equation (5). Since the intersecting areas result in the reduction of the coating volume fraction, the coating layers in the RVE should be thickened to compensate the areal loss due to the interferences. After the fiber distribution process is completed, the thicknesses of all the coatings are uniformly increased by a factor of , defined as the square root of the ratio between the coating area in the RVE and the reference coating area, . is the measured coating volume fraction, and is the side length of the RVE. As the coatings expands, the intersecting areas sequentially grow, which may lead to an additional loss in the entire coating area. Therefore, this compensation process should be iteratively conducted until the area ratio, , becomes one, which indicates the equality between the total coating area considering the intersections and the reference coating area. Figure 5 shows the overall process for the construction of an RVE model with the randomly distributed coated fibers.
When determining the size of an RVE model, it is known that the model size should be sufficiently smaller than the characteristic length of composites and sufficiently larger than that of constituents to obtain the credible estimation of effective composite properties [37]. As discussed by Kanit et al. [38], the RVE size determines the precision of the property predictions. In other words, when using small RVEs, one needs to generate a sufficiently large number of RVEs to obtain reliable computational results. The main purpose of this study does not focus on the precise estimation of homogenized elastic properties but focuses on the fracture behaviors and postpeak responses, which are greatly affected by the spatial arrangement of fibers. Thus, the size of RVE in this study is set to contain 20 fibers, which is comparable with the previous researchers [2, 34, 35] who have also studied microstructural fracture behaviors with nonuniform fiber arrangements. Figure 6 shows the RVE models of the same size with the same fiber and coating volume fractions, but different fiber distributions. They will be utilized to study the effect of the fiber distribution on the microstructural behavior of the CMCs.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4. Microstructural Fracture Behavior of CMCs
4.1. Finite Element Model
Finite element analysis (FEA) using a commercial software package, ABAQUS/Explicit, is performed to study the effect of the random fiber geometries on the microstructural fracture behavior of the CMC material based on the RVE models shown in Figure 6. All the RVEs have 20 coated fibers in the same size squares with the side length of 90.65 μm. The dimensions and locations of the coated fibers are different between the eight RVEs, but the fiber and coating volume fractions of all the models are identical. The volume fractions of the fibers and coatings are 0.35 and 0.06, respectively, as reported in Section 2. The RVE models are meshed with over 18,000 three-node plane strain triangle elements (CPE3) as shown in Figure 7. The corresponding mesh size is sufficiently fine such that any further refinement does not significantly improve the accuracy of FEA results. Transverse tensile loading and pure shear loading as shown in Figure 8 are applied to the RVE models to study the effect of the randomness in the fiber geometries on the mechanical responses of the CMC materials. Periodic boundary conditions (PBCs) are utilized to implement the loading conditions.
(a)
(b)
4.2. Periodic Boundary Conditions
The PBCs are kinematic constraints enforced on nodal degrees of freedom between two paring nodes on opposite sides. A set of equations for the application of the PBCs to a three-dimensional (3D) model can be found in [39]. For the two-dimensional (2D) RVE models as shown in Figure 9, the constraint equations for the two vertical edges are expressed asand the remaining equations for the horizontal edges are
The DOFs on the nodes at the vertices are doubly constrained by using both equations (21a) and (22b). By eliminating common terms in the redundant constraints, the PBCs for the DOFs at the vertices are reduced towhere and are the DOFs of the nodes on the edges in the and directions, respectively. , , and can be considered as macroscopic strains that determines the global deformation of the 2D structure. As can be seen in Figure 9, the global strains or the global deformation fields, and , are assigned on the two master nodes 1 and 2 and the DOFs on the edges are associated with them through equations (21a)–(23d). Table 1 shows the displacement setup for the tensile and shear loading conditions considered in the present study. The periodic boundary conditions expressed in equations (21a)–(23d) are implemented into the finite element model (FEM) through “” option of ABAQUS.
4.3. Smeared Crack Approach
The microstructural fracture behavior of the CMC material is studied using the smeared crack approach (SCA) proposed by Rots et al. [40]. SCA does not explicitly model a discontinuity due to a crack but incorporates its effect into the framework of continuum mechanics by degrading stiffness in a postpeak softening regime. The effect of a crack is incorporated by assuming that the total strain increment, , is decomposed into the continuum strain, , and the cracking strain, :
When a crack is initiated in an element at an angle of as shown in Figure 10, the strain increment due to the crack in the local coordinate system associated with the crack normal and tangent vectors can be transformed into the global strain increment such thatwhere is the transformation matrix between the local and global coordinate systems and is the crack strain increment under the local coordinate system. It should be noted here that the present study considers Mode I crack only by assuming that the crack is opening due to the maximum principal stress [41]. When a crack is considered in the principal plane with the maximum principal stress, the mixed mode fracture behavior is inherently taken into account.
The local stress increment, , can also be related to the global stress increment, , through the transformation matrix:
The constitutive relation between and is assumed to bewhere the secant stiffness matrix, , can be identified from the traction separation law depending on the loading status. For example, when the fractured element is being unloaded, the secant stiffness, , as shown in Figure 11 is entered into equation (27).
When the continuum stiffness of the element in Figure 10 is , the constitutive relation between the continuum strain in equation (24) and the continuum strain increment , can be expressed as
Using equations (25)–(27), equation (28) can be written in terms of the crack stiffness matrix and the transformation matrix it reads:or in a conventional form of a damage parameter model:where is the identity tensor and can be considered as a damage index that degrades the continuum stiffness . Note that contains the crack stiffness matrix determined from the current crack situation as illustrated in Figure 11, and thus the effect of the crack is incorporated into the continuum model. When the maximum principal stress of the element exceeds the cohesive strength, a crack is initiated in the element and the elemental stiffness is being degraded through . Before the crack strain increases to the final deformation, , as defined in Figure 11, the element is considered to be being damaged due to the crack growth. When the crack strain reaches , the crack is considered to be fully opened and the element loses its stiffness completely.
In presentation, SCA used for coating the layer and the matrix material only. The fracture behavior of the fiber is not considered here since the fiber is a highly fracture-resistant material compared to the other two constituents. The fiber is simply modeled as an isotropic linearly elastic material. The prepeak behavior of the coating and the matrix is also assumed to behave like an isotropic linearly elastic material. The material and fracture properties used in the present numerical study are listed in Table 2. The SCA is implemented into the present FE model through a user subroutine option of ABAQUS/Explicit.
4.4. Effective Elastic Properties
Before the present modeling approach is applied to study the microstructural fracture behavior of the CMC materials, its validation is first quickly assessed through the predictive capability of effective elastic properties. Table 3 compares the transverse elastic stiffnesses obtained from an experiment and numerical models. Table 4 lists the constituent properties used by the numerical models. Here, we generate four RVEs with different fiber arrangements to predict the transverse property. Totally 20 fibers of the same diameter and interphase thickness are used for all the four RVE models. The average and the corresponding standard deviation of the four results from the present model are reported in Table 3. The present modeling approach shows the best agreement with the experimental data in terms of the relative error. When the relative error is defined in Table 3, the experimental value is used as a true value. The average value of the present model is used to compute the error in Table 3.
4.5. Transverse Tensile Loading Case
Figure 12(a) shows a typical stress-strain response of an RVE under the transverse tensile loading condition as described in Table 1. The crack area fraction, defined as the ratio of the total crack area to the square RVE area, is also displayed in Figure 12 as a function of the axial strain. Note that the dramatic changes of the crack area are closely related to the load drops in the stress-strain curve as shown in Figure 12(a). The progressive damage extents associated with the stress-strain curve are also illustrated in Figures 12(b)–12(e). The stress-strain response is initially linear but starts to exhibit plastic-like behavior until the first peak at stage A. The first load drop is caused by the fully opened crack in the matrix material as shown in Figure 12(b). In addition to the fully matured crack, there are multiple growing cracks in the coating and material phases at stage A. Indeed, the multiple existing cracks, including the fully opened one before it is fully opened, degrade the stiffness of the entire RVE structure, and this appears as the “pseudoplastic” behavior in the global stress-strain curve before the first peak. The second load drop at stage B results from the multiple cracks found on the right-hand side of the RVE model as illustrated in Figure 12(c). The first fracture dissipates highly accumulated elastic energy around the region, and the stresses in the entire RVE area are redistributed. While the stress relaxation occurs around the area of the first fracture, the other side is already very highly stressed at stage B. As a result, multiple cracks emerge almost simultaneously on the opposite side of the first crack as shown in Figure 12(c). As the multiple cracks grow, they join together to form a longer transverse crack, which eventually splits the RVE into two pieces at stage D.
(a)
(b)
(c)
(d)
(e)
Figure 13 shows the stress-strain responses of the RVE models shown in Figure 6. The corresponding final fracture patterns are displayed in Figure 14. As can be seen in Figure 13, the initial slopes of all the curves are the same since the different RVE models have the same fiber volume fraction. Beyond the elastic regime, the stress-strain curves start diverging and the postpeak responses are distinctively different from each other. It is also interesting to note that, for all the RVEs, the linear-to-nonlinear transitions occur slightly above 50 MPa, which is close to the cohesive strength of the coating. This indicates the phase that fractures first and responsible for the plastic-like stress-strain response.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.6. Pure Shear Tests
The shear response of the RVE models are also studied using the shear loading condition described in Table 1. Shear strains of up to 0.8% is applied to all the RVE models in Figure 6. Figure 15(a) shows a typical shear stress-shear strain response of an RVE together with the evolution of a fractured area. Every sudden increase in the area, mainly due to the emergence of a new crack and/or the growth of existing cracks, is again directly associated with the load drops of the stress-strain curve. Fracture patterns at which the RVE exhibits particular stress-strain behavior are also depicted in Figures 15(b)–15(e). In contrast to the case of the transverse tensile loading, the matrix is damaged over almost the entire area since the shear load is mainly carried by the low-stiffness material. Note that undamaged regions are colored blue, green indicates a crack initiation, and red is used to designate the region where a crack is fully opened.
(a)
(b)
(c)
(d)
(e)
The first load drop at stage A in Figure 15 is caused by the initiation of a crack, of which the orientation is perpendicular to the loading-direction diagonal of the square model as shown in Figure 15(b). It is reiterated here that, in the present study, only Mode I fracture is considered through the SCA. It is assumed that an element is considered to be fractured when the maximum principal stress of the element exceeds the corresponding cohesive strength and the crack in the element is oriented perpendicular to the principal direction. For the pure shear loading condition considered here, it is obvious that the principal directions of each element in the matrix phase are closely aligned with the loading-direction diagonal of the RVE model. Consequently, the first crack shown in Figure 15(b) is oriented perpendicular to the particular direction. Between the stages A and B, more cracks are initiated while the existing cracks have grown. At stage C in Figure 15(d), multiple short cracks have emerged in addition to the lengthy cracks, resulting in the severe fluctuation in the shear stress-shear strain curve after the first load drop. The orientations of the major cracks found at stage C are still closely associated with the loading-directional diagonal of the RVE model. A significant load drop is observed between C and D. The existing cracks grow and join together even with a newly initiated long crack as shown in Figure 15(e). It is interesting to note that the new crack is developed in the direction perpendicular to the existing ones. This implies that the principal directions of the elements in the RVE model have been altered significantly since stresses are redistributed whenever a crack is present.
The shear stress-shear strain responses of all the RVE models in Figure 6 are shown in Figure 16. The final fracture patterns associated with each stress-strain curve are displayed in Figure 17. Similarly, as with the results of tensile tests, when the shear strain is less than 0.001, the initial stiffness of the RVEs are almost identical mainly due to the same fiber volume fraction. After the shear strain exceeds 0.001, the stress-strain curves start showing a slight deviation from each other, indicating that the materials are being damaged. The shear responses of the RVE models significantly differ from each other beyond the shear strain of 0.003. The fracture patterns at the final shear strain value in Figure 16 are also very complicated and unique depending on the random fiber distributions in the RVE models as shown in Figure 17. It should be noted here that the cracks are diagonally oriented for the case of the shear loading.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
5. Conclusions
A micromechanical model for the transverse cross section of CMCs has been developed to study the effect of nonuniform fiber distributions on the fracture behavior when the model is subjected to tensile and shear loading, respectively. A kinematic algorithm creates an RVE model having nonuniformly distributed coated fibers, while the corresponding geometric parameters and volume fractions are controlled to be the same as the measured data. The virtual tests based on FEA are performed, and the RVEs showed fairly the same responses in the linear region. However, a great variety is observed in the nonlinear and postpeak regions mainly due to the different rates of damage propagation and cracking patterns affected by the fiber arrangement.
Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Korean Government (Ministry of Trade, Industry and Energy, MOTIE, and Defense Acquisition Program Administration, DAPA) through Institute of Civil-Military Technology Cooperation.