The transverse strength of fiber-reinforced composites is a matrix-dominated property whose accurate prediction is crucial to designing and optimizing efficient, lightweight structures. State-of-the-art analytical models for composite strength predictions do not account for fiber distribution, orientation, and curing- induced residual stress that greatly influence damage initiation and failure propagation at the microscale. This work presents a novel methodology to develop an analytical solution for transverse composite strength based on computational micromechanics that enables the modeling of stress concentration due to Repre- sentative Volume Elements (RVE) morphology and residual stress. Finite element simulations are used to model statistical samples of composite microstructures, generate stress-strain curves, and correlate statistical descriptors of the microscale to stress concentration factors to predict transverse strength as a function of fiber volume fraction. Tensile tests of thin plies validated this approach for carbon- and glass-reinforced com- posites showing promise to obtain a generalized analytical model for transverse composite strength prediction. Keywords – Transverse tensile strength, analytical modeling, computational micromechanics, damage mechanics 2 1 Introduction Designing advanced composite parts for engineering and structural applications is extremely time- consuming as it often involves expensive and lengthy experimental campaigns or complex and time-intensive computational modeling [1]. It is paramount to expedite the composite design and certification process to meet their ever-increasing demand. Analytical models can be used as efficient predictive tools to estimate composite response from constituent material properties, especially in problems involving many design vari- ables for optimization [2, 3], if their relative error with respect to experimental measurements is contained. State-of-the-art models employ either physics-based formulations or empirical relations to correlate the overall composite response to their constituent structures and properties. Analytical models provide significant ad- vantage of simple implementation and computation over experiments and high-fidelity computational models for structural design and optimization [4]. The literature presents numerous analytical models that can estimate the elastic property of polymer matrix composites (PMCs) from their constituent properties. For instance, classical theories such as rule- of-mixture (ROM), modified rule-of-mixture (MROM) and enhanced models such as the continuous periodic fiber model (CPFM), the concentric cylindrical assembly (CCA), Mori-Tanaka (MT), and the Bridging Model (BM), have been widely employed to predict composite elastic response [3–14]. However, there are only a few micromechanical strength theories that can accurately predict the composite strength from its constituent properties [2, 6, 8, 15–17]. Huang et al. [6] developed the BM to estimate the stress field inside a single-fiber repeating unit cell (RUC) with square packing. They introduced a stress concentration factor (SCF) within the RUC to compute the in-situ matrix strength which was then used in conjunction with explicit micromechanical formulae to evaluate the transverse composite strength. Their model presented an average error of 40% when compared with transverse strengths of nine unidirectional (UD) composite laminates. Vignoli et al. [2] developed an elasticity-based model for transverse strength prediction. They established analytical solutions to estimate the stress distribution around a single-fiber in an infinite medium. Combining them with finite element (FE) approach, adjustments functions were defined to determine the stress distribution in a more realistic square-packed RUC of several fiber volume fractions. Finally, the transverse strengths were predicted based on the matrix cavitation model. Their predictions presented an average error between 25% and 40% when compared to experimental data as well as other simplified analytical models. These models, among others, were developed based on a simplified mathematical representation of a composite microstructure, or a single-fiber square RUC [2, 3, 8, 18, 19]. Thus, they failed to account for inter-fiber interactions and the resulting stress concentration distribution that is observed in a multi-fiber representative volume element (RVE). Furthermore, the models did not consider the effect of process-induced in-situ matrix property variation and residual stress generation on the transverse composite strength. It is well known that random fiber packing and residual stress generation not only influence the failure initiation in composite microstructures, but also drive failure propagation under continued loading [4,20–27]. Existing analytical models ignore these effects and 3 therefore present a lower correlation with transverse composite strength tests. Computational micromechanics- based simulations can fill these knowledge gaps and address the fundamental challenges associated with analytical models in order to enhance their prediction accuracy. Over the past few decades, advances in high-performance computing, finite element methods (FEM) and novel physics-based constitutive modeling have facilitated the design and development of PMCs through computational methods [28–30]. Computational micromechanics is an emerging field that leverages these advances to model and analyze RVEs of composite microstructures by subjecting them to various thermo- mechanical boundary conditions [30–50]. Such analysis can provide valuable insights into the fundamental mechanisms that influence the composite response at the constituent-level, thereby overcoming the limitations associated with experimental quantification [30]. By explicitly modeling the constituent fibers and matrix, phenomena such as stress concentration due to random fiber packing and fiber-to-fiber proximity, that induce matrix cracking under transverse tensile loading conditions and significantly influence the transverse compos- ite strength, can be studied with relative ease and efficiency. FE-based process modeling simulations, which are informed by accurate and comprehensive material characterization and employ phenomenological and con- stitutive relations, can facilitate the quantification of residual stress build-up within composite RVEs during manufacturing. Virtually curing and loading RVEs in transverse tension can provide insight on the impact of process-induced residual stresses and random fiber distribution on transverse composite strength. The knowl- edge gained from computational micromechanical analysis of composite RVEs can be leveraged to establish a processing-microstructure-property relationship and address the aforementioned fundamental challenges to develop an accurate computational micromechanics-based analytical model for transverse strength. The objective of this work is to present a novel methodology to develop a computational micromechanics- based analytical model for transverse tensile strength prediction. In this study, FE-based, computational micromechanical models are generated to analyze a wide array of composite microstructures for their transverse composite strengths. It is hypothesized that random fiber packing and the resulting stress concentration induce failure within the microstructures, which affects their transverse strength. Based on this hypothesis, a generalized expression for transverse strength is discussed. Critical parameters that influence the transverse composite strength are highlighted. The concept of a stress concentration factor induced by a single-fiber RUC is presented. The influence of relative fiber arrangement (distance between fiber centers and their orientation) on the SCF is then quantified through FE simulations of two-fiber model. Leveraging the insights from the two-fiber SCF study, a statistical approach to quantify the SCF in multi-fiber microstructures is presented. Finally, a material system-specific analytical model for transverse composite strength is developed. It is shown that insights from computational models can inform the development of mathematical relations to correlate the microscale physical mechanisms such as stress concentration to composite failure and predict transverse strength of multi-fiber composite microstructures. This manuscript is organized as follows: the computational micromechanical modeling approach is detailed 4 in Section 2, the development of the analytical model is described in Section 3 and the concluding remarks are discussed in Section 4. 2 Computational Micromechanics Computational micromechanics requires efficient yet sufficiently large models to capture the physical mechanisms that induce failure. In this study, converged RVEs including 50 fibers were analyzed based on the results of a statistical size-effect study performed previously by Shah et al. [32]. The modeling details and virtual manufacturing and testing procedures are discussed in the following sections. 2.1 Modeling of composite microstructures In this study, two sets of RVEs comprising a random dispersion of IM7 carbon (fiber diameter, df= 6µm) and E-glass fibers ( df= 14 µm) in an epoxy matrix, respectively were generated using a random RVE generator developed by Stapleton et al. [51]. A commercial epoxy system, EPIKOTETMResin MGS RIMR 135 with EPIKURETMCuring Agent MGS RIMH 1366 (henceforth referred to as RIM R135-H1366), was chosen as the matrix material due to its comprehensive thermo-mechanical material property dataset reported in [31]. Several composite microstructures, corresponding to various fiber volume fractions (0 .25≤vf≤0.75) were considered, as illustrated in Figure 1 and summarized in Table 1. To improve the prediction accuracy and to account for statistical variations due to random fiber architecture, five distinct replicates with random fiber distribution were analyzed for each value of vf. Perfect bonding was assumed between the fibers and the matrix. Due to their nature, the constituent IM7 carbon fibers were modeled as transversely-isotropic while the E-glass fibers were modeled as isotropic solids. The thermo-mechanical properties for both constituent fiber materials are summarized in Table 2. The matrix material was modeled as isotropic. The evolution of the matrix thermo-mechanical properties with the degree of cure ϕand the processing temperature Tduring virtual manufacturing was defined by experimentally-determined material-specific models described in [31]. These properties, for a fully cured matrix ( ϕ= 1), are also listed in Table 2. Failure was admissible only in the matrix material when the maximum principal stresses exceeded its critical strength. 5 (a) IM7 carbon fibers RVEs (b) E-glass fiber RVEs Figure 1: Various realizations of (a) IM7 carbon fibers (CF) and (b) E-glass fibers (EGF) embedded in RIM R135-H1366 epoxy matrix with varying fiber volume fractions 25% ≤vf≤75%. 6 Table 1: Details of composite microstructures analyzed in this study. Material RVE ID RVE length No. of fibers vf (µm) (-) CF-25 75.2 50 0.25 CF-35 63.5 50 0.35 CF-45 56.0 50 0.45 IM7 carbon fiber CF-55 50.7 50 0.55 (df= 6µm) CF-60 48.5 50 0.60 CF-65 46.6 50 0.65 CF-75 43.4 50 0.75 EGF-25 175.5 50 0.25 EGF-35 148.3 50 0.35 E-glass fiber EGF-45 130.8 50 0.45 (df= 14 µm) EGF-55 118.3 50 0.55 EGF-65 108.8 50 0.65 EGF-75 101.3 50 0.75 Table 2: Thermo-mechanical properties of the constituent fibers and epoxy matrix. Property Units IM7 carbon E-glass RIM R135-H1366 Density, ρ [kg/m3] 1780 2550 1200 Axial Modulus, E11 [MPa] 276000 73000 2482 Transverse Modulus, E22=E33 [MPa] 19500 73000 2482 In-plane Poisson’s ratio, ν12=ν13 [-] 0.28 0.22 0.37 Out-of-plane Poisson’s ratio, ν23 [-] 0.25 0.22 0.37 In-plane Shear Modulus, G12=G13 [MPa] 70000 30000 905.8 Out-of-plane Shear Modulus, G23 [MPa] 7800 30000 905.8 Critical Tensile Strength, σcr [MPa] – – 64.1 Mode I Fracture Energy, GIC [J/m2] – – 0.001 Axial CTE, α11 [K−1] -0.54 ×10−65×10−673.9×10−6 Transverse CTE, α22=α33 [K−1] 10.08 ×10−65×10−673.9×10−6 Thermal Conductivity, k [W/m-K] 5.4 1.2 0.245 Specific Heat, cp [J/kg-K] 879 800 1600 The virtual manufacturing and mechanical testing procedures were carried out with the commercial FE code Abaqus/STANDARD, supplemented with user-written UMATHT and UMAT subroutines, the imple- mentations of which are detailed in [31,32] and summarized in the following section. Each RVE was meshed with C3D8T elements (eight-node, brick elements with temperature degrees of freedom). Flat boundary con- ditions (FBCs) were chosen for both virtual manufacturing and mechanical testing, as illustrated in Figures 2a and 2b, respectively [32]. 7 (a) Virtual manufacturing (b) Virtual mechanical testing (c) Cure and thermo-mechanical property evolution for a prescribed cure cycle Figure 2: Illustration of flat boundary conditions (FBCs) prescribed to the RVEs during: (a) virtual manu- facturing and (b) virtual mechanical testing analyses; (c) evolution of the degree of cure ϕ, glass transition temperature Tgand the elastic modulus Efor a prescribed cure cycle. 2.2 Virtual manufacturing The virtual manufacturing analysis replicated the complete manufacturing of composites by accounting for (a) the kinetic progression of the cure which was quantified by ϕand (b) the evolution of the in-situ thermo-mechanical matrix properties which led to residual stress generation. The progression of cure, for a prescribed temperature profile, was defined by the Kamal-Sourour kinetic model [52], dϕ dt= (K1+K2ϕm) (1−ϕ)n(1a) Ki=Aiexp∆Ei RT∗ i = 1,2 (1b) 8 where K1andK2are Arrhenius rate functions; m= 0.4 and n= 1.5 are dimensionless, real modeling parameters; A1= 3.6×109sec-1andA2= 0.01245 sec-1are pre-exponential factors; ∆ E1= 85.3 kJ/mol and ∆ E2= 11.1 kJ/mol are activation energies; T∗is the absolute temperature in degrees Kelvin; Ris the universal gas constant. The kinetic constants were experimentally determined as reported in [31]. The kinetic model was solved simultaneously with the 3-dimensional Fourier heat transfer model in a coupled-temperature displacement analysis, carried out in Abaqus/STANDARD with user-subroutine UMATHT, to predict the progression of cure and the temperature distribution resulting from the exothermic heat of reaction during curing. The prescribed cure cycle and the computed degree of cure is presented in Figure 2c. The evolution of the in-situ matrix properties and residual stress generation was modeled with UMAT user subroutine. For a given cure state, the instantaneous material properties were computed using experimentally- determined material-specific models reported in [31]. For instance, the variation of the in-situ matrix elastic modulus Ewith ϕandT(see Figure 2c) was defined by E=A2(ϕ) +A1(ϕ)−A2(ϕ) 1 + exp η−η0(ϕ) dη(ϕ) (2) where the dimensionless phenomenological parameters η,η0, dη,A1andA2were determined as a function of ϕandTby the following expressions η= 1−T∗ T∗r(3a) η0=−0.860ϕ+ 0.955 (3b) dη= 0.070ϕ+ 0.041 (3c) A1(ϕ) =−10−13exp(34 .55ϕ) (3d) A2(ϕ) = 1471 .3ϕ+ 801 .39 (3e) where T∗ ris the reference temperature in degrees Kelvin. Similarly, the evolution of the matrix glass transition temperature Tϕ g(see Figure 2c), post-gelation chemical shrinkage ϵϕ shand post-gelation coefficient of thermal expansion αϕ iwere determined using the following equations Tϕ g=T0 g+λϕ 1−(1−λ)ϕ T1 g−T0 g (4a) ϵϕ sh=  0, ϕ ≤ϕgel β(ϕ−ϕgel), ϕ > ϕ gel(4b) αϕ i=  7.93×10−5, T ≤Tϕ g [a−b(ϕ−ϕgel)]×10−5, T > Tϕ g(4c) 9 where T0 g=−53.2◦C and T1 g= 85.6◦C are glass transition temperatures of uncured and fully cured matrix; β= 0.111 is the shrinkage coefficient; ϕgel= 0.78 is the gel point (see Figure 2c); λ= 0.2373, a= 21.21 and b= 18.25 are dimensionless fitting parameters. An instantaneous linear-elastic constitutive model was then employed to model the residual stress generation in the material during curing σi(t) = [Cij(t)] ϵtotal j(t)− ϵth j(t) +ϵsh j(t) δj where,  δj= 1 if j = 1, 2, 3 δj= 0 if j >3(5) where i and j are Voigt notation indices; ϵtotal j(t),ϵth j(t) and ϵsh j(t) are the total, thermal and chemical shrinkage strains, respectively; Cij(t) is the stiffness matrix as a function of the time of cure; σi(t) is the accumulated residual stress governed by the development of the in-situ matrix elastic modulus and the chemical and thermal strains experienced by the matrix. In this study, the curing matrix was prescribed a constant strength and fracture toughness as listed in Table 2. As the matrix material, subjected to the temperature profile shown in Figure 2c, cured, its elastic mod- ulus developed as illustrated in Figure 2c. This, combined with the matrix chemical shrinkage and thermal mismatch between the constituent fiber and matrix, resulted in self-equilibrating residual stress generation. The contour plot of the residual stresses (maximum principal), at the end of the cure, is shown for one rep- resentative IM7 carbon and E-glass fiber microstructure in Figures 3a and 3b, respectively. These RVEs, which consisted of a random dispersion of 50 fibers with a volume fraction of 0.6 and 0.55, respectively registered a maximum residual stress generation of 21.48 MPa and 24.58 MPa. The same microstructures reported a volume-averaged residual stress of 6.38 MPa and 7.49 MPa, respectively. These maximum and volume-averaged residual stress values differed considerably between RVEs owing to the different constituent fiber material properties, their unique fiber packing topology and the varying fiber volume fractions. Between the IM7 carbon fiber RVEs with vf= 0.6, there was a maximum and volume-averaged end-of-cure residual stress of 20 .93±0.46 MPa and 6 .406±0.073 MPa, respectively. The E-glass fiber RVEs, with vf= 0.55 reported a maximum and volume-averaged end-of-cure residual stress of 24 .79±1.1 MPa and 7 .89±0.44 MPa, respectively. 10 (a) IM7 carbon fiber RVE ( vf= 0.60) (b) E-glass fiber RVE ( vf= 0.55) (c) Stress versus strain plots Figure 3: Contour plots showing the end-of-cure residual stresses (maximum principal) from virtual man- ufacturing analysis, maximum principal stress before the onset of damage from virtual mechanical testing and damage field at the end of the virtual mechanical testing in representative (a) IM7 carbon fiber RVE (vf= 0.60) and (b) E-glass fiber RVE ( vf= 0.55); (c) plot of stress versus strain from virtual mechanical testing of IM7 carbon and E-glass fiber RVEs. 2.3 Virtual mechanical testing The processed RVEs at the end of the previous step were subjected to transverse tensile mechanical loads by prescribing a displacement boundary condition as illustrated in Figure 2b. The objective here was to compute the transverse composite response of the RVEs, namely the transverse tensile strength σ+ 22as a 11 function of their fiber volume fraction. Failure in the matrix was modeled with a previously developed progressive damage model [32], based on the crackband theory [53]. The maximum principal stress criterion was utilized to determine failure initiation in the matrix. A traction-separation law, governed by the fracture energy, was employed to define the post- peak softening behavior of the damaging material once the critical fracture stress was reached. The critical mode I energy release rate GICwas given by GIC=hηZ¯ϵη f 0¯ση 11¯ϵη 11d¯ϵ (6) where ¯ ση 11and ¯ϵη 11are the maximum principal stress and strain values in element η, respectively; ¯ ϵη fis the value of ¯ ϵη 11which corresponds to a zero stress state on the post-peak stress versus strain plot; hηis the characteristic length of the element ηthat preserves mesh objectivity by prescribing a normalized value of GICfor each element, such that gη IC=GIC/hη. A scalar damage factor Dηwas computed to degrade the damaged element compliance components using Dη= 1−σcr Em(¯ϵη f−¯ϵη init)¯ϵη f ¯ϵη 11−1 (7) where ¯ ϵη initis the value of ¯ ϵη 11when the initiation criterion (¯ ση 11≥σcr) is satisfied; Emis the undamaged Young’s modulus of the matrix. The damage parameter could take values between zero and one, where Dη= 0 meant no damage had occurred. By contrast, a maximum damage level of one corresponded to a zero-stress state on the post-peak stress versus strain plot. Also, healing was inadmissible. Once the damage factor was computed, the relevant components of the compliance matrix were degraded [32,54]. The progressive damage formulation was modeled in Abaqus/STANDARD solver with user-written subroutine UMAT. The matrix strength σm cr and a scaled-down fracture energy GIC, corresponding to the sub-micron length scale, were prescribed to the material as listed in Table 2. The stress versus strain plots of the five IM7 carbon fibers RVEs ( vf= 0.6) and E-glass fiber RVEs (vf= 0.55), subjected to post-manufacturing transverse loads, are presented in Figure 3c. In this figure, the RVEs exhibited an initial linear-elastic response which was followed by a drop in the stress. This drop was associated with damage initiation in the matrix when the local stresses exceeded the matrix strength σm cr. The contour plot of the maximum principal stress before the onset of damage in the representative microstructure is shown in Figure 3a for IM7 carbon and Figure 3b for E-glass fiber RVEs, respectively. Local regions with high stress concentration (corresponding to warmer colors), where failure initiation was anticipated, were clearly visible in these plots. The contour plot of the end-of-loading damage variable D, which was computed using Equation 7 and used to degrade the elements stiffness, is presented in Figures 3a and 3b, respectively for IM7 carbon and E-glass fiber RVEs. It was evident from these figures, that stress localized in regions of dense fiber packing which led to failure initiation and local microcracking. These cracks eventually coalesced 12 into a large crack that propagated through the microstructure resulting in a two-piece failure. The peak stress in the stress versus strain plots was regarded as the transverse composite strength σ+ 22. For the representative IM7 carbon and E-glass fiber microstructures shown in Figure 3, σ+ 22= 24.03 MPa and σ+ 22= 17.49 MPa, respectively. The analysis was repeated for all RVEs illustrated in Figure 1 and their replicates. The transverse com- posite strength predictions from these analyses are summarized in Figure 4 as a function of their fiber volume fraction for IM7 carbon and E-glass fiber RVEs. Each data point in these plots was obtained by averaging the numerical predictions from five distinct replicates of an RVE of a given fiber volume fraction. The relevant standard deviations from the average values were represented by the corresponding error bars. It was clear from these figures, the fiber volume fraction manifested a strong influence on the transverse strength of the composite. The transverse composite strength predictions for IM7 carbon fiber RVEs increased steadily with the fiber volume fraction from σ+ 22= 21.23±1.37 MPa for vf= 0.25 to σ+ 22= 28.65±0.71 MPa for vf= 0.75. Similar trends were observed for E-glass fiber RVEs where the transverse strength predictions increased from σ+ 22= 15 .94±2.61 MPa for vf= 0.25 to σ+ 22= 20 .74±0.32 MPa for vf= 0.75. Furthermore, the fiber packing manifested a strong influence on the transverse strength. Random and densely packed fibers induce local stress concentration in composites leading to premature failure initiation and local microcracking, which typically manifests itself as large scatter of bulk composite strength. This was evident from the large standard deviations in the predicted strength values for each volume fraction in Figure 4. In particular, RVEs with lower fiber volume fraction ( vf= 0.25,0.35) manifested a relatively larger scatter associated with their pre- dicted strength values. This suggested that the converged RVE size of 50 fibers, which was defined based on a fiber volume fraction vf= 0.55, may not be sufficiently large and representative of composite microstructures with lower fiber volume fractions. To define an appropriate RVE size for such cases, a separate statistical size effect study is warranted, which can be computationally expensive due to the larger RVE sizes that need to be analyzed. Notwithstanding that, the standard deviations for these cases were found to be within 10% of the mean values, which is within the accepted limits for numerical and physical testing. Furthermore, the predicted values were directly compared with the respective experimental strength values obtained from in-house transverse testing of thin IM7 carbon and E-glass fiber composite laminates. The predicted strength values were in excellent agreement with the experimental values as presented in Figure 4. 13 (a) IM7 carbon fiber RVEs ( df= 6µm) (b) E-glass fiber RVEs ( df= 14 µm) Figure 4: Summary of the transverse composite strength predictions from the virtual testing of (a) IM7 carbon fiber and (b) E-glass fiber RVEs composite microstructures as a function of their fiber volume fraction. 3 Development of the Analytical Model Modeling [21,25,32] and testing [55] of composites have shown that matrix failure initiates in regions of high stress concentration as a result of dense fiber packing. This concentration of stress not only depends on the constituent fiber and matrix properties but also the random fiber distribution, residual stress generation and the global fiber volume fraction that influences the transverse composite strength. This work proposed an analytical model that directly accounted for these parameters in Equation 8 σ+ 22=σ+ 22 Em, Ef 22, νm, νf, σm cr,SCF+ 22 (8) specifically, the transverse composite strength σ+ 22was expressed as a function of the constituent fiber/matrix elastic properties Em,Ef 22,νm,νf, the in-situ matrix strength σm crand a stress concentration factor SCF+ 22. Results from Section 2 also showed that the relative distance between the neighboring fibers and their ori- entation with respect to the loading direction significantly influenced stress localization between them which eventually resulted in failure initiation. Therefore, the SCF was divided into two contributions, one that accounted for the effect of the fiber distribution (distance between the fiber centers and their orientation), SCF I, and the second that accounted for the effect of the fiber volume fraction, SCF II.SCF+ 22was defined by, SCF+ 22=SCF I(δ, θ)×SCF II(vf) (9) where, δis the shortest distance between a fiber pair while θis the orientation of that shortest distance with respect to the loading direction and vfis the fiber volume fraction. The composite transverse strength can be 14 calculated starting from the definition of stress concentration as shown in Equation 10 σ+ 22=σm cr SCF I×SCF II(10) Solving for σ+ 22required the determination of the SCFs. A novel computational approach to calculate the SCF IandSCF IIfor a multi-fiber RVE is proposed below. 3.1 Stress concentration factor A stress concentration factor is defined as the ratio of the stresses in the proximity of a discontinuity within a continuum and the far-field/applied stress averaged over a boundary surface. In case of composite microstructures, a fiber inclusion within a matrix is a discontinuity that acted as a stress-riser. A simple linear-elastic analysis of a single IM7 carbon fiber embedded in an infinite RIM R135-H1366 medium yielded an SCF = 1.45. This meant the stresses in the vicinity of a fiber were 1.45 times higher that the applied stresses. A similar value of SCF = 1.54 was obtained for E-glass fibers embedded in the same RIM R135-H1366 epoxy matrix. While such results from a single-fiber RUC have been extensively used to develop analytical models for transverse composite strength, they fail to provide useful information about the non-uniform stress field in a multi-fiber RUC. In such cases, the neighboring fibers influence the stress distribution around the fiber being considered, which then affects the stress concentration surrounding it. This evaluation of a global SCF in a multi-fiber microstructures becomes increasingly challenging since the non-uniform stress field is not only influenced by several neighboring fiber interactions within a certain proximity δbut also their relative orientation with respect to the loading direction θ. To determine the SCF in multi-fiber microstructures, a unique stress concentration analysis was proposed in this study. First, the influence of introducing a neighboring fiber on the local SCF was investigated with a two-fiber model (see Section 3.2). Following this, a correlation was established between the geometrical variations in a multi-fiber microstructure and the SCF by virtue of the nearest neighbor statistical descriptor (see Section 3.3). 3.2 Influence of fiber arrangement on SCF The influence of fiber arrangement on the SCF was investigated through FE simulations. Several two- fiber models, as illustrated in Figures 5 and 6, were generated in Abaqus/STANDARD with 0◦≤θ≤90◦ and 6 .1µm≤δ≤26µm for IM7 carbon fiber models while 0◦≤θ≤90◦and 14 .1µm≤δ≤34µm for E-glass fiber models, respectively. Assuming a linear elastic material response, these models were subjected to transverse mechanical loads by prescribing a displacement boundary condition. The SCF in each case was determined using the classic definition, SCF I=σmax σ∞(11) 15 where, σmaxis the maximum principal point stress and σ∞is the far-field/applied principal stress. The computed SCF Ifor various combinations of δandθare shown in Figures 5a and 5b for IM7 carbon fibers and Figures 6a and 6b for E-glass fibers, respectively. For a given value of δ= 7.2µm in Figure 5a, the highest stress concentration was observed when the fiber centers were oriented parallel to the loading direction ( θ= 0◦). TheSCF Irapidly decreased from 2.52 to 1.45 as θincreased from 0◦to 90◦. By contrast, the SCF Igradually reduced from 3.13 to 1.56 in Figure 5b as the distance between the fiber centers δincreased from 6.2 µm to 15µm while θ= 0◦. A consistent trend was observed for E-glass fiber models as illustrated in Figures 6a and 6b. The influence of relative fiber arrangement on the local SCF was clearly evident from these figures. By analyzing an exhaustive combination of δandθ, an interpolation function for SCF I(δ, θ) was developed (see Figures 5c and 6c), SCF I(δ, θ) =a11+a12θ+a21δ+a22θδ+a23δ2(12) where, a11−a23are fitting parameters specific to the material systems considered in this study. These values for both IM7 carbon and E-glass fiber models are summarized in Table 3. Note, these values were different for the two fiber materials analyzed in this study and that their values may further vary for other material systems. This function, defined by Equation 12 and the corresponding fitting parameters, were utilized in the subsequent section to determine the SCF Iof each multi-fiber microstructure analyzed in the previous section as a function of their geometrical variations. 16 (a) Variation in SCF Iwith θwhen δ= 7.2µm (b) Variation in SCF Iwith δwhen θ= 0◦ (c) Variation in SCF Iwith θandδ Figure 5: Variation in SCF Iin IM7 carbon fiber microstructures with (a) fiber orientation θand (b) distance between the fiber centers δ(fibers not shown); (c) interpolation function expressing SCF Ias a function of θ andδ. 17 (a) Variation in SCF Iwith θwhen δ= 15.2µm (b) Variation in SCF Iwith δwhen θ= 0◦ (c) Variation in SCF Iwith θandδ Figure 6: Variation in SCF Iin E-glass fiber microstructures with (a) fiber orientation θand (b) distance between the fiber centers δ(fibers not shown); (c) interpolation function expressing SCF Ias a function of θ andδ. 18 Table 3: Material system-specific fitting parameters used for the analytical model predictions. Parameter [-] IM7 carbon E-glass SCF Single fiber RUC 1.45 1.54 a11 Equation 123.698 13.3 a12 -1.872 ×10-2-7.035 ×10-2 a21 -2.123 ×10-1-7.86×10-1 a22 -9.742 ×10-42.508×10-3 a23 -4.645 ×10-31.278×10-2 b Equation 132.99×10-11.336×10-1 c 5.916 8.449 d -5.882 -8.811 3.3 SCF in multi-fiber microstructures Several statistical descriptors have been reported in the literature to quantify the geometrical variations within the microstructures such as local fiber volume fraction distribution using Voronoi tessellation, cluster analysis with Delaunay triangulation, second-order intensity function [26, 27, 56–62]. Since the goal of this study is to quantify the local geometrical variations in the multi-fiber microstructures analyzed in Section 2 and to establish a direct correlation with the global SCF I, the nearest neighbor statistical descriptor was utilized which is based on the short-range fiber interactions. The definition of this descriptor is illustrated schematically in Figure 7. For a given fiber center iin the microstructure, the algorithm identified its closest neighbor and measured the distance between the fiber centers δialong with its relative orientation θiwith respect to the loading direction. This process was repeated for all fiber centers in the microstructure. The mean shortest distance ¯δand the corresponding mean orientation ¯θwere computed for each microstructure. The global SCF Ifor each of these microstructures was then computed with the help of the interpolation function SCF I¯δ,¯θ (see Equation 12), where the mean shortest distance ¯δand orientation ¯θvalues from the statistical analysis were used. These values for each microstructure analyzed in this study are summarized in the appendix Tables A1 and A2. 19 Figure 7: Schematic showing the working of the nearest neighbor distribution statistical descriptor. By substituting into Equation 10 , the neat matrix strength σm cr= 64 .1 MPa (see Table 2), composite strength predictions σ+ 22of microstructures from Section 2 and the SCF Icomputed using Equation 12, SCF II was determined for each microstructure analyzed previously (see appendix Tables A1 and A2). The variation ofSCF IIfor each microstructure as a function of its fiber volume fraction is presented in Figure 8 for IM7 carbon and E-glass fiber models. Note that each data point on the plot is an average of five microstructure renditions with random fiber packing. A non-linear least-square regression fitting of this variation yielded an expression for SCF II, SCF II=b+c√vf+dvf (13) where b−dare fitting parameters specific to the material systems and microstructures considered in this study. These values for both IM7 carbon and E-glass fiber models are summarized in Table 3. Equations 12 and 13 were substituted in Equation 10 to obtain the final expression for the transverse composite strength. The analytical model, thus developed, accurately predicted the transverse strength of IM7 carbon and E-glass fiber composite microstructures for a wide range of fiber volume fractions and random fiber distributions. For instance, in Figure 9a, the model predicted close to the neat matrix strength of σm cr= 64.1 MPa where vf= 0 which is representative of neat matrix. Subsequently, a sudden drop in strength predictions was observed with the introduction of fibers due to stress concentration. This drop was followed by a steady increase in the composite strength with the fiber volume fraction. A good correlation was established between the numerical model predictions, experimental measurements and the analytical model predictions as demonstrated in Figure 9a. Consistent trends were obtained form E-glass fiber/RIM R135-H1366 microstructures as seen in Figure 9b. 20 (a) IM7 carbon fiber RVEs ( df= 6µm) (b) E-glass fiber RVEs ( df= 14 µm) Figure 8: Variation in SCF IIwith the fiber volume fraction of (a) IM7 carbon fiber and (b) E-glass fiber RVEs. (a) IM7 carbon fiber RVEs ( df= 6µm) (b) E-glass fiber RVEs ( df= 14 µm) Figure 9: Comparison of the analytical model predictions for transverse composite strength with the numerical model predictions and experimental data for (a) IM7 carbon and (b) E-glass fiber RVEs. It should be noted that Equations 10, 12 and 13 are material system-specific and cannot be employed to predict the transverse composite response of a different constituent material system. This limitation is associated with the underlying assumption that the influence of the constituent fiber/matrix elastic properties on the SCF was implicitly accounted for in the data fitting process. As summarized in Table 3, a different set of fitting parameters were obtained for E-glass fiber microstructures. Thus, to better quantify the influence of constituent material properties on the SCF and to develop a generalized analytical mode, independent of any fitting parameters, explicit correlations between the constituent properties and the SCF must be defined. Additionally, the effect of process-induced residual stresses on the transverse strength predictions was implicitly accounted for during the data fitting process where composite strength of the virtually manufactured microstructures were used to inform the analytical model development. To explicitly determine this effect 21 and to account for it in the improved analytical model, further investigation into the development of residual stress and the influence of random fiber distribution is warranted. It has been shown that fiber can influence how stress is generated during processing and redistributed during subsequent loading which drives failure initiation and propagation in composites. A good understanding of these physical phenomena is necessary before its influence on the transverse composite strength can be explicitly modeled. Once its influence is quantified, a dedicated parameters can be introduced to capture the effect of process-induced residual stress on the transverse composite strength prediction. Finally, the present model was developed using the nearest neighbor distribution (based on short-range interaction) statistical descriptor to quantify the random fiber distribution and determine failure initiation in composite microstructures. However, it was observed that post- initiation, failure propagation through the microstructure not only depended on the fiber closeness near the initiation site but also on the local fiber distribution (local fiber volume fraction) that influenced residual stress generation during manufacturing and stress redistribution during subsequent mechanical loading. To quantify such geometrical variations and to fully capture their effect, more sophisticated statistical descriptors must be assessed. Despite these limitations, the proof of concept to develop a computational micromechanics-based analytical model for transverse composite strength prediction is promising and warrants further investigation. 4 Conclusions This work presented a novel methodology to develop an analytical model for transverse composite strength for carbon and glass reinforced epoxy resin using a numerical and statistical approach. Computational micromechanical analysis, involving virtual manufacturing and testing, were carried out on composite mi- crostructures with varying fiber volume fractions to predict their transverse composite strength. A novel stress concentration factor for multi-fiber microstructures were defined based on the nearest neighbor statisti- cal descriptor. A material system-specific analytical model was developed for transverse strength predictions based on the computational modeling results. It was shown that computational micromechanics can provide insight on the physical mechanisms that influenced the transverse composite strength, which are otherwise extremely challenging to quantify experimentally. Accounting for such physical mechanisms, a mathemati- cal relation was established, which linked microscale failure in composites to these physical mechanisms, to effectively predict the transverse composite strength. Acknowledgments This material is based upon the work supported partially by the Air Force Office of Scientific Research and the National Science Foundation under grant number IIP – 1826232. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. The authors would like to thank NASA for their support of this 22 research under grant numbers 80NSSC19K1246 and NNX17AJ32G. No potential competing interest was reported by the authors.