1. Introduction
Automation of the manufacturing process of composites may well be the next important cost efficiency gain in industrial applications (Beukers et al. [
1]). Mass produced textile fiber performs that can be molded into near net shapes will lead to composite structures that have longer fatigue lives. There is an urgent need for robust and high fidelity computational tools that can replace costly testing associated with certification of composites for aerospace structural applications. Driven by the necessity of certification by analysis, integrated computational materials engineering (ICME) has emerged as a vibrant research activity in the past few years. Robust, mechanics based models that possess the correct physics at each stage of a production cycle will be needed to ensure maturity of ICME based computational frameworks. Summary of the major challenges related to ICME of composites can be found in Arnold et al. [
2].
We focus on an important aspect of the ICME of fiber reinforced composites pertaining to composite unit cell at the micro-scale. Following earlier work (Heinrich et al. [
3]), we address the issue of unit cell size that is needed for assessing the stiffness and strength of such unit cells when subjected to mechanical loads. The study includes virtual modeling that includes the development of damage and cracking as the unit cells are loaded. During processing of thermoset fiber reinforced polymer composites (FRPCs), the matrix undergoes shrinkage due to network formation and thermal strains due to thermal mismatch with surrounding fibers, which give rise to internal stresses. The state of the matrix at the end of curing can be altered by the presence of fibers and also by the form of the cure cycle. Song et al. [
4] showed that the use of bulk matrix properties (i.e., measured “virgin” matrix properties) in numerical predictions of compression response of 2D triaxially braided composite representative volume element can lead to erroneous results. They reported that the computed compressive strength using the virgin matrix properties were higher than experimentally measured strength. Such a discrepancy in strength prediction using the so called “virgin” matrix properties arises because matrix that cures in the presence of the fibers is noticeably different from the virgin matrix. Such a matrix is called in situ matrix, which has the effective properties of the matrix that takes into account imperfections caused in the matrix during the curing process. The post-cure in situ matrix properties in terms of the matrix shear stress-shear strain response, can be extracted via an inverse analysis through uniaxial tensile loading of a ±45
laminate, which is often convenient in engineering analysis of processed composites (Ng et al. [
5]). The in situ shear stress-shear strain response can also be obtained from a torsion test. Thus, good knowledge of the state of matrix in the curing FRPC structure is important for ICME to be successful for FRPCs. Moreover, for a particular FRPC system, the optimal cure cycle must be identified such that the cured structure has the desired strength and stiffness. In the present investigation, we use micromechanics to study the influence of unit cells of certain size, which contain randomly packed fibers. The unit cells are chosen to be of different size (preserving the fiber volume fraction but containing different numbers of randomly packed fibers). The effects of curing of unit cells of different sizes is subject of ongoing work—that is, the presence of residual stresses in the composite due to processing is neglected in the present analysis.
Figure 1 qualitatively shows the behavior of a composite specimen at the global scale concurrently with the behavior of fiber-matrix constituents during initiation and propagation of damage up to subsequent failure. Although global damage is observed macroscopically, which manifests as matrix cracking resulting in a reduction of tangent stiffness, the matrix may have failed completely locally. Due to the inability to model the entire fiber-matrix microstructure at the component (global) level, the material properties at the global scale are usually the homogenized properties obtained from the constituents, using the properties of stiffness and strength from the unit cell. In a multiscale analysis setting using finite elements, a simple unit cell representation is often used to represent the effective behavior of the global structure at the finite element’s integration points. Therefore, it is important to first understand how the choice of the unit cell size affects its overall stiffness and strength representation of much larger volumes of fiber-matrix composite. Thus, the goal of this work is to show whether or not stiffness and strength of unit cells are controlled by its size.
In this study, the matrix tensile strength used is similar to what is representative of macroscopic coupon level measurements, while the Mode I fracture toughness
has been lowered by approximately two orders of magnitude compared to macroscopically measured values (for example, from a double cantilever beam test or a compact tension test with specimens in the tens of mm scale). The use of macroscopically observed values of
in the microscale simulations encountered in this study would lead to physically unrealistic response. Recently, Qiao and Salviato [
6] proposed a two-scale cohesive traction-separation law. Based on measurements, they showed that the fracture strength of thermoset polymers at the microscale is noticeably higher than the macroscopically observed fracture strengths, whereas the fracture toughness,
, associated with Mode I cracking is up to two orders of magnitude lower than what is observed in macroscopic specimens. Their observation has also been supported by coarse-grained molecular dynamics (MD) studies on the tensile behavior of thermoset polymers ( Yang and Qu [
7]). This MD study showed that sharp reduction in cohesive strength occurs due to cavity initiation and coalescence phenomenon at the sub-micron scale associated with the crack-opening displacements. The subsequent macro-scale reduction in cohesive strength and toughness is due to re-alignment of the polymer chains in relatively less densely packed volumes. This suggests that the fracture toughness and strength measured at the macroscopic scale are not appropriate for use in microscale fracture simulations.
2. Modeling
In reality, the microstructure of fiber-reinforced composites exhibits randomness in packing (
Figure 2), which in turn influences the mechanical properties (Bulsara et al. [
8], Moulinec and Suquet [
9]). In contrast, the baseline or the idealized case is the one where fibers are arranged in a hexagonally packed manner. One of the questions we can ask is that for a given fiber volume fraction, how much do the transverse fracture properties (strength and toughness) of the randomly packed configuration deviate from the perfect arrangement. Another aspect is determining the number of fibers (or unit cell size) for given (fixed) fiber volume fraction sufficient enough to represent the bulk transverse fracture properties. The problem of finding the smallest unit cell that will represent bulk fracture properties is affected more by the geometrical details of the constituents compared to corresponding problem for stiffness, where the bulk properties are governed by the average properties of the constituents.
In an earlier work by Heinrich et al. [
3] that addressed the influence of packing on the residual stress buildup during cure, the randomness in packing was characterized by Voronoi tessellation (
Figure 3). The deviation from hexagonal packing was described in terms of the area of the Voronoi cells, the centers of which corresponded to the centers of fibers placed in the unit cell. The present investigation will follow a similar approach of characterizing fiber packing. For each unit cell size and for a given fiber volume fraction, several fiber packing renditions will be created.
Our prior work (D’Mello et al. [
10]) on transverse loading of virtually cured randomly packed unit cells indicates that packing has a strong influence on transverse tensile strength (
Figure 4). However, in this investigation, apart from influence of unit cell size, we also aim to find the influence of randomness in fiber packing and matrix fracture properties under transverse tension. Under tension, the crack band model described in the next section, will be used to model damage/failure propagation in the unit cells. The output from these simulations will provide transverse tension fracture properties—that is, unit cell strength and fracture toughness. Next step in the analysis is finding correlation between unit cell size (for a given fiber volume fraction) and the fracture properties and quantifying randomness associated with fracture properties. As shown in
Figure 4, taken from D’Mello et al. [
10], the effect of curing has a marked influenced on the unit cell strengths. This behavior is due to the presence of inherent self-equilibrating residual stresses in the unit cells on account of curing. When such unit cells are loaded in the transverse direction, the overall strain required for the matrix elements to enter failure is much lower. Moreover, the global stress recorded is also lower due to account of pre-existing tensile residual stress in the matrix.
In the present work, random fiber unit cells with three different sizes are considered: unit cells having 6, 24, 40 and 180 fibers respectively. Each of these unit cells will be loaded in the transverse direction (two-direction) under tension. Commercially available finite element software ABAQUS/Explicit complemented by user subroutines are used for these numerical simulations. Each of the unit cells are modeled in 3D, using linear hexahedral C3D8 elements. The fiber diameter is taken to be eight microns and the dimension of the unit cell in the fiber-direction is kept fixed at 0.2 times the fiber diameter. The mesh density is taken such that the circumference of the fiber is modeled with approximately 38 C3D8 elements. For the unit cells with six fibers, three different fiber volume fractions (
) of 0.4, 0.5 and 0.6 are considered, to study the variation of fiber volume fraction on the stiffness and strength, whereas for the unit cells having 24, 40 and 180 fibers, a fiber volume fraction of 0.5 is considered in order to study the effect of size on the unit cell stiffness and strength. For each type of unit cell, 10 different realizations of randomly packed fiber configurations are considered to quantify the scatter in the mechanical response under transverse tension. Each type of unit cell used in the present work along with the coordinate frame indicated are shown in
Figure 5. Carbon fibers are modeled as transversely isotropic solids with properties
235,000 MPa,
15,000 MPa,
MPa,
= 27,000 MPa,
and
. The matrix has Young’s modulus
MPa and Poisson’s ratio
. During tensile loading in the transverse direction (two-direction), the unit cells are subjected to flat boundary conditions, wherein the edges of the unit cell are allowed to expand or contract under loading but are constrained to remain flat throughout the loading process.
A crack band model based on the one proposed by Bazant and Oh [
11] is used to model failure in the matrix. Methodology of crack band model for unit cells have been addressed in our earlier paper (D’Mello and Waas [
12]) but will be repeated here for convenience. This model assumes that once the critical fracture stress
is reached locally in the matrix, microcracks are formed and this effect is smeared over an element. The maximum principal stress criterion is used to determine the failure initiation in the matrix. However, other failure criteria can also be used for transverse compression, based on experimental observations and data that is relevant for modeling transverse compressive failure. For example, Pineda et al. [
13] have used Mohr–Coulomb failure criterion for compressive-shear failure in fiber-reinforced composite unit cells using finite element analysis as well as by using high-fidelity generalized method of cells (Aboudi et al. [
14]). In the present study for unit cells under transverse tension, when the maximum principal tensile stress in the matrix exceeds the critical value
, the traction–separation law controls the behavior of the damaging material as shown in
Figure 6 and the stiffness of the matrix is reduced and consequently tracked using the secant modulus. When a traction–separation law controls the behavior of the damaging solid, the response is that of a non-continuum in the softening regime. Under mode I cracking, the energy dissipated during failure is the critical mode I energy release rate
, given by
where,
and
is the maximum principal stress and principal strain, respectively, the maximum separation
where
corresponds to the strain corresponding to complete failure of the material (accompanied by complete loss of stiffness). Here,
h is the characteristic element length that preserves mesh objectivity (see Jirasek and Bazant [
15]), by prescribing a normalized value of
for each element such that
. Consequently, the value of
equals the area under the
law as shown in
Figure 6. In the finite element procedure, the finite element sizes were kept more or less uniform in all the unit cells considered. The characteristic length
h is computed for each element via Abaqus’s internal characteristic element length CELENT variable (Abaqus [
16]). The value of
was chosen to be 0.002 N/mm in all the computations, whereas the tensile fracture strength
was chosen to be 55 MPa. From the crack band model formulation, the damage parameter due to softening is denoted as
D where
. For a material with initial matrix stiffness
now in the softening region of the traction-separation law,
D is computed as:
where
is the current maximum principal strain. Thus, the case
corresponds to no damage in the finite element, the case
corresponds to damage but no two-piece failure in the element, while the case
indicates failure with complete loss of stiffness. Thus, the quantity
is a measure of stiffness reduction in the damaging matrix with as a fraction of the original stiffness
. Note that although we are dealing with a 3D solid, the stress–strain equation for the traction separation law is provided in the 1D form because the relevant stress and strains are assumed to be arising from an equivalent definition of stress and strain. In this case, the equivalent stress and strains are taken as maximum principal values as described in Bazant and Oh [
11]. The damage parameter
D is calculated based on this assumption and is used to degrade the 3D solid.
The notion of an representative volume elements (RVEs) is most appropriate for finding stiffness (elastic) properties of the bulk. However, once failure has set in (i.e., in the regime past the peak in the stress–strain response of the unit cell), the concept of an RVE in the traditional sense is lost. Failure in realistic bulk materials occurs in a localized manner. Thus, we can no longer assume that failure process is replicated in all periodic arrays of material in the bulk. The existence and determination of RVE has been addressed in detail by Gitman et al. [
17], who demonstrated that an unit cell cannot be found for a material in the softening regime. Later, Nguyen et al. [
18] have also addressed the issue of existence of the unit cells for softening quasi-brittle materials with random microstructure.
It is also instructive to point out and differentiate the quantities related to various sizes used in the finite element simulations relevant to the discussion of the crack band model. First is the typical length of the elements that must be maintained in order to achieve stress convergence, especially in zones where stress gradients are high, such as in the vicinity of fibers. In the present study, is proportional to circumferential length of the fibers, that is , where is the fiber diameter. This length is independent of the actual length scale of the structure but prescribed in proportion to fiber geometry, so that the computed stresses around fibers have converged. Next, independent of , there are additional lengths and which must ensure that the tangent strain–softening modulus, for modes I and II cracking respectively, associated with the crack band model is not negative. These two lengths depend on the length scale of the element. Here, for mode I cracking and for mode II cracking, where G is the shear modulus and is the critical shear strength for mode II cracking. For use in finite elements, recommended limits for and are half of what is prescribed by the expressions. Although stresses may have converged with the use of converged , we must ensure that and , so that the size of elements are suitable for use in the crack band model.