Multiscale Modeling of Fiber Fragmentation Process in Aligned ZnO Nanowires Enhanced Single Fiber Composites

A three-dimensional multiscale modeling framework is developed to analyze the failure procedure of radially aligned zinc oxide (ZnO) enhanced single fiber composites (SFC) under tensile loading to understand the interfacial improvement between the fiber and the matrix. The model introduces four levels in the computational domain. The nanoscale analysis calculates the size-dependent material properties of ZnO nanowires. The interaction between ZnO nanowires and the matrix is simulated using a properly designed representative volume element at the microscale. At the mesoscale, the interface between the carbon fiber and the surrounding area is modeled using the cohesive zone approach. A combination of ABAQUS Finite element software and the failure criteria modeled in UMAT user subroutine is implemented to simulate the single fiber fragmentation test (SFFT) at the macroscale. The numerical results indicate that the interfacial shear strength of SFC can be improved up to 99% after growing ZnO nanowires on the fiber. The effect of ZnO nanowires geometries on the interfacial shear strength of the enhanced SFC is also investigated. Experimental ZnO nanowires enhanced SFFTs are performed on the fabricated samples to validate the results of the developed multiscale model. A good agreement between the numerical and the experimental results was observed.


Multiscale Modeling of enhanced Sfcs with Aligned Zno nanowires
The multiscale analysis and modeling of the enhanced SFCs are described in this section, which is divided into four length scales: (i) single ZnO nanowire at nanoscale, (ii) ZnO reinforced epoxy unit cell at microscale, (iii) interfacial enhancement of aligned ZnO arrays on the carbon fiber at mesoscale, and (iv) the SFFT and fiber fractures at macroscale, as shown in Fig. 2. The modeling results at each length scale are implemented as input for the next length scale. nanoscale analysis of single Zno nanowire. A core-shell structure is used to model ZnO nanowires and to evaluate the material properties of ZnO nanowires with different geometries. According to the theoretical approach proposed by Chen et al. 18 and verified with experiments, the size-dependent Young's modulus of ZnO nanowires with diameters in the range of 17-450 nm can be calculated using Eq. 1.   where E 0 is the modulus of the bulk ZnO material. Considering a core-shell structure of the nanowires as a composite wire, E s is the surface modulus, r s is the depth of the shell and D is the nanowire diameter. The obtained value of E 0 is 139 GPa, r s is 4.4 nm, and E s /E 0 is equal to 1.50 18 . Besides, various typical lengths of nanowires in the range of 2.5-14 µm are taken into account in this study, which can control the enhanced ZnO/epoxy coating layer as shown in Fig. 1.
Microscale modeling of Zno reinforced epoxy unit cell. In order to study the interfacial enhancement of hybrid SFC, the coating layer (radially aligned ZnO) is modeled as nanowires reinforced epoxy as shown in Fig. 2b,c. A small section of the ZnO/epoxy coating layer with the periodic unit cell can demonstrate the properties of the whole geometry. According to Fig. 2b, a representative volume element (RVE) can be considered in which, a single ZnO nanowire is embedded in the matrix. The maximum ZnO volume fraction in the coating layer can be achieved with the compact radial alignment of the nanowires orthogonal to the carbon fiber as demonstrated in Fig. 2c. Considering an element consists of a single ZnO with the length of (L ZnO ) reinforced the epoxy shown in Fig. 2b, the maximum volume fraction of ZnO can be defined as Eq. 2. The detail of this calculation is described in the Appendix.
where d f is the diameter of the carbon fiber and L ZnO is the length of the ZnO nanowire. The effective material properties of the coating layer can be estimated by the homogenization of this RVE based on the theory of continuum micro-mechanics. The bonding between ZnO and epoxy in this RVE is assumed to be perfect. According to the theory of mechanics, the stress-strain equations for an anisotropic material is defined by Eq. 3.
where E ij is Young's modulus, G ij is the shear modulus, and ν ij is the Poisson's ratio of the orthotropic material in the ij = 1, 2, 3 directions, respectively. Assuming 1 is the longitudinal direction, 2 is the in-plane transverse direction, and 3 is the out of plane transverse direction. Considering the volume of the RVE as V RVE , the average strain ( ij ε ) and stress (σ ij ) of the composite can be calculated from Eqs. 5, 6.
In this homogeneous RVE, the total energy of the system (U) is obtained from Eq. 7, while the stored energy in the system (U * ) can be found from Eq. 8.
According to the equilibrium state of the strain energy in the system, U * − U = 0. By applying the different boundary conditions and displacements in the RVE model, and using equilibrium principle of energy, the matrix www.nature.com/scientificreports www.nature.com/scientificreports/ of elastic compliance can be calculated 20 . More information about extracting the effective material properties of a square cross-section RVE using the continuum mechanics can be found in the ref. 21 . The approach provided by Omairey, et al. 22 is used in ABAQUS for obtaining the effective material properties of the ZnO nanowires reinforced epoxy. The simulated model of this analysis is shown in Fig. 3a considering the perpendicular orientation of the carbon fiber and ZnO depicted in Fig. 3b, the obtained material properties should be updated based on the coordinate system of the carbon fiber (upper letter) and the ZnO nanowire (lower letter).
Mesoscale analysis of the fiber/matrix interface. The adhesion bonding between the carbon fiber and the surrounded area is modeled at mesoscale based on the cohesive zone model and the cohesive elements. In this method, the interface is modeled as a very thin layer. The bilinear traction-separation law is considered to model the stress-displacement of the interface, which has been widely used in the FEM modeling of the carbon fiber reinforced epoxy 23,24 . According to this model, when the stress in the interface layer reaches its strength, the progressive damage initiates and the stress degrades linearly to model the separation of the interface 25 as shown in Fig. 4a. The maximum stress in the interface is defined as nominal or shear interface strength (σ max ) and the fracture energy (G Ic ) is described as the area under the bi-linear curve. The stress-displacement relation for different parts of the bilinear traction-separation can be defined by Eq. 9 26 .
where D is the damage parameter and K is the stiffness in the interface. The D parameter is defined with Eq. 10 using the maximum stress criteria for modeling the damage in the interface.  www.nature.com/scientificreports www.nature.com/scientificreports/ The values of the bilinear parameters can be calculated with the experiments. In the case of carbon fiber and epoxy matrix, 40-60 MPa has been reported for the σ max 24,27 , while the value of G Ic has been obtained between 100-230 J/m 2 28 . The values of the maximum stress and the fracture energy is considered as (σ max = 50 MPa, G Ic = 100 J/m 2 ) in this study. Calculation of a specific value for the stiffness with the experimental or analytical analysis is a challenge, and a wide range can be found in the literature. Accordingly, the mathematical approach is implemented here for calculating the interface stiffness based on the material properties of the fiber and the surrounded materials. Different radial variation of the interface modulus across its thickness has been explored in the literature with the boundary condition defined as Eq. 11 29 .
where E f and r i show Young's modulus and the radius of the fiber, while E m and r i describe the modulus of the matrix and the radius of the interface, respectively. The effective stiffness of the interface (E i ) can be defined as the average of the modulus along with the thickness written in the form of Eq. 12 30,31 . Based on the power-law variation shown in Fig. 4b, the interface modulus is written as E i = Ar B , and the average interface modulus in this model can be written as Eq. 13 29,30 .
Macroscale analysis of Sfft. The adhesion bonding between fiber and matrix can be evaluated via the SFFT. The applied tensile load is transferred to the fiber through the interface between matrix and fiber by the shear transferring mechanism. By increasing the load, the fiber breaks into the small fragments after the stress in the fiber exceeds its strength. Each fragment can carry the load to form independent segments and breaks again. This phenomenon continues to form a progressive fiber failure to the point that the segments become too small to break described as saturation level. According to the shear stress model introduced by Kelly-Tyson 32 , the constant interfacial shear strength can be calculated using Eq. 14.
where r f and σ f * are the radius and the ultimate strength of the fiber, respectively. The parameter of (L c /2) is the critical distance from the fiber segment end in which the stress reaches the fiber strength. According to the model provided by Ohsawa et al. 33 , the critical length (L c ) is defined as a function of the average length of the segments (L i ) in the saturation state with n fragments as shown in Eq. 15. The mechanism of the stress on different fragments of the fiber is depicted in Fig. 5.
In order to simulate the SFFT, the failure of the fiber is simulated using the Tsai-Wu failure criterion 34 . The original form of Tsai-Wu expression for orthotropic materials can be written in the form of Eq. 16 35 . www.nature.com/scientificreports www.nature.com/scientificreports/ Assuming the fiber as transversely isotropic material with the tensile strength of σ t * , compressive strength of σ c * , and shear strengths of τ * with the same values in the longitudinal and transverse direction, the F ij coefficient can be written as Eq. 17. The magnitude of F is called as a reference to the failure situation in which, the material is safe when F has a magnitude of smaller than 1, and the failure occurs if F equals to 1. The F factor is defined and tracked using UMAT user subroutine programmed in FORTRAN. When the damage occurs, the related element cannot carry the load while other elements in the model are active, which breaks the fiber in different segments. The general steps of the damage analysis are depicted in Fig. 6.

feM Analysis of Sfft
A 3-D model of single carbon fiber embedded in the epoxy is simulated. In order to minimize the computational cost, a quarter of a 3-D composite beam is modeled considering two symmetric planes. Based on the symmetry on the x-plane, the displacement in the x-axis (U x ), rotation around y-axis (UR y ), and rotation around the z-axis (Ur z ) are set to zero. The displacement in the y-axis (U y ), rotation around the x-axis (UR x ), and rotation around the z-axis (Ur z ) are zero for the y-symmetric plane. The axial displacement load is applied to the x-faces as shown in Fig. 7. A thin layer of interface is considered between the fiber and matrix based on the cohesive zone concept to model the adhesive bonding. The FEM analysis of the SFC performed in the ABAQUS package consists of two stages, including bare fiber and enhanced fiber. In the bare SFC, carbon fiber with 7 μm diameter is placed in the center of a cubic matrix with a width of 0.14 mm and a length of 0.5 mm. It has been shown that the interface thickness has a negligible effect on the utilization of the cohesive elements while a value close to zero is suggested 36,37 . Hence, an interface thickness of 0.01 μm is connected to the fiber and epoxy using "tie constraint". In the enhanced fiber model, the carbon fiber with the same geometry is coated with a layer of ZnO/epoxy composite embedded in the same geometry of matrix to form the four geometry phases. The ZnO nanowires are aligned radially on the surface of the fiber. Accordingly, the thickness of this coating layer is controlled by the length of the ZnO. The material properties of the carbon fiber and epoxy used in this analysis are described in Table 1 38 . The tensile, compressive, and shear strength in this table are shown with σ t * , σ c * , and τ * , respectively. The simulated SFC with the defined boundary conditions is depicted in Fig. 7. An 8-node linear brick element with improved surface stress visualization (C3D8S) is utilized for the matrix, ZnO/epoxy layer, and the fiber. An 8-node 3-D cohesive element (COH3D8) is used to model the interface.

Multiscale Modeling Results of Sfft
Homogenization of the ZnO/epoxy coating layer. The appropriate RVE is considered to calculate the effective material properties of the ZnO/epoxy composite layer. The maximum volume fraction of 57.87% is calculated from Eq. 2 for the ZnO with a length of 2.5 μm. The cubic RVE based on this volume fraction is simulated in ABAQUS (Fig. 3a) considering different diameters and the FEM homogenization is performed. The modulus of the related nanowire diameters extracted from Eq. 1 is implemented in the homogenization analysis. The results for nanowires with diameters of 17, 75, 150, 450 nm are depicted in Table 2. Besides, assuming a constant diameter of 17 nm, the effective material properties of four different nanowire lengths (L = 2.5, 5, 10, 14 μm) were calculated as shown in Table 3. Due to the radial arrangement of nanowires on the fiber, the volume fraction of ZnO in the ZnO/epoxy layer is reducing by increasing the length of the nanowires according to Eq. 2. Based on the relation between the ZnO and the carbon fiber longitudinal axis described in Fig. 3b, the appropriate material properties calculated from the RVE homogenization is used to simulate the ZnO/epoxy coating layer in the enhanced SFC.
Stress distribution along the fiber. Considering the material properties and the cohesive zone model described in the previous sections, the stress transferred into the fiber from the matrix is calculated in different situations. In order to evaluate the mechanism of the load transferring, the stress distribution along the fiber   www.nature.com/scientificreports www.nature.com/scientificreports/ length is evaluated for the case of bare SFC. This trend can be compared with the shear-lag theory proposed by Cox 39 . According to this approach, the stress distribution on the fiber can be extracted using Eqs. 18, 19.
where r f is the radius of the fiber and ε m is the applied strain and V f is the volume fraction of the fiber. Besides, E f and E m are Young's modulus of the fiber and the matrix, respectively. The results calculated by the theory and the FEM for the 0.5% applied strain is shown in Fig. 8. It can be observed that the stress reaches its maximum value at the center of the fiber and degrades to the zero value at the two ends. Besides, the maximum value of the stress in the fiber calculated by FEM is fairly correlated with the theory.
fiber fragmentation results. The mesh dependency analysis was first performed to find the proper mesh size for the FEM of SFFT. Hence, the maximum axial stress on the fiber under 0.6% applied strain is obtained for different mesh density as shown in Table 4. Since the relative difference in the FEM results for the 18836 total element numbers is 0.6%, this mesh density is chosen for the analysis. The maximum axial stress on the fiber is increased by improving the applied displacement load on the composite and the failure occurs when the damage model described in the previous section is satisfied. After the failure, the fiber breaks into two different segments each of which can carry the load separately up to the critical state. The stress distribution along the fiber axis at different state of fracture, including before the first break, after the first break, after the third break, and after the five breaks is calculated for the enhanced SFC with d ZnO = 17 nm and L ZnO = 2.5 μm as shown in Fig. 9. It can be observed that the fracture occurs in the middle of each segment, while the stress distribution in the smaller segments shows similar trend along the fiber axis.
The constitutive failure is continuing up to the state that the segment length in the fiber is small for the stress to reaches the fiber strength. In order to explore the effect of the enhanced interface to the load-carrying properties of the fiber, the density of the fiber fragmentation in the bare fiber is compared to the results of the fiber with the ZnO (d = 17 nm and L = 2.5 μm) coating layer as illustrated in Fig. 10a. According to this figure, the number of fractures in the enhanced fiber (n = 19) is almost doubled of the bare fiber (n = 9), which means more load is transferred to the fiber through the enhanced interface. Besides, according to this improvement, the first fiber fracture occurred at the lower applied strain in the case of enhanced fiber compared with the bare one. The total number of breaks at the saturated state is used for calculation of the interfacial strength from Eqs. 14, 15. Accordingly, the interfacial shear strength for the bare carbon fiber is 95.5 MPa compared to 190 MPa for the  Table 4. Mesh density effects on the fiber's maximum stress.
enhanced fiber. In other words, by adding the ZnO/epoxy composite layer on the surface of the carbon fiber, the interfacial shear strength is enhanced by 99%. In addition, the stress-strain curve for the SFC sample is shown in Fig. 10b for the bare and the enhanced fiber. The degradation in the slope of the stress-strain curve indicates the fracture in the fiber. The area where the first fracture happened is marked as an example for both samples. Besides, according to the stronger interface, the first fracture occurs at lower strain compared to the bare sample.
Effect of the ZnO nanowire geometries. The effect of ZnO nanowires geometry on the enhancement of the interface between fiber and matrix was evaluated. In this regard, the effective material properties for the nanowire with various diameter in the range of 17-450 nm extracted from the homogenization section was implemented in this analysis. The density of the fiber fragmentation calculated from FEM for the different diameter is illustrated in Fig. 11. The interfacial shear strength for each diameter was calculated as shown in this figure. It can be seen that the thin ZnO nanowires result in an increased number of fragmentations (stronger interface). The nonlinear trend indicates the large impact on the interface enhancement when thinner nanowire is applied. The interfacial shear strength for the 17 nm diameter is 190 MPa compared with the 114 MPa for the nanowires with 300 nm diameter. In other words, the interface between fiber and epoxy is stronger when smaller ZnO is used. In addition, the length effect of the nanowires on the interface was explored. The effective material properties of the coating layer with ZnO length of 2.5, 4, 10, and 14 μm is implemented for FEM based on Table 3. By extracting the number of fragmentations, the interfacial shear strength is calculated for each geometry as it is shown in Fig. 12a. It can be observed that the interface is weakened almost linearly by increasing the length of the nanowire. This behavior can be also observed in Fig. 12b in which the strain at the first fracture is shown for different lengths. From this figure, the first fracture occurred earlier in the shorter nanowires based on the stronger interface that leads to transferring more load from the matrix to the fiber.

experiments
SFFT experiments were carried out in this study to validate the interfacial properties of SFC using either bare carbon fiber or ZnO nanowires enhanced carbon in SFC. Atomic layer deposition (ALD) and hydrothermal methods were employed to synthesize the aligned ZnO nanoarrays on carbon fiber fabrics. The detailed nanowire synthesis procedure has been reported in the author's recent publication 40 . The ZnO nanowires have an average length of    Fig. 13a. The mold cavities were then filled with an epoxy resin (EPON 862/Curing Agent 9553: 100/17.6 by weight) and cured at the elevated temperature for 6 hours at 120 °C. The specimens were tested in a screw-driven micro-tensile machine (Fig. 13b) under a constant displacement rate of 0.2 mm/min, and the fragmentation process was recorded using an optical microscope and polarized light. The specimen was unloaded after the saturation point was reached, and the number of fragments was measured by taking pictures under polarized light. Figure 14(a) shows a typical optical micrograph after tensile testing taken under polarized light. The experiments show that growing the ZnO nanowires has a negligible effect on the fiber's tensile strength. Fiber breaks, matrix cracks, and interfacial debonding can be observed from the micrograph. The epoxy matrix is optically isotropic but becomes anisotropic due to local stress concentration. In addition, the interfacial debonding between fiber and matrix are highlighted in bright light due to the optically anisotropic matrix in the debonding area. In the FEM, the strain redistribution around the fiber fracture area results in stress/strain concentration on the matrix. These critical areas obtained from the nine fragmentations of the FEM model can be compared with the damage zone of the fiber/matrix in the microscope image as shown in Fig. 14a,b. The experimental results of the number of fractures with the applied strain are shown in Fig. 14c. It can be observed that the average number of fragmentations are increased in the ZnO nanowires enhanced SFC. This result demonstrates that growing ZnO nanowires on the fiber results in a strong interface enhancing the bonding between the carbon fiber and the epoxy matrix. The stronger interface leads to the more efficient load transferring from the matrix to the fiber in the SFFT for the enhanced fiber. Hence, the surface area of each fiber segment has enough shear load to transfer to the fiber to exceeds the strength and causes more fractures.
The average number of fragmentations is 33 and 20 for the hybrid and bare fiber, respectively. The experiments show 88% increase in the interfacial shear strength of the enhanced fiber compared to the bare fiber assuming the same fiber tensile strength, and growth of 0.98 µm in the diameter of the enhanced fiber. The numerical analysis of the radially aligned ZnO nanowires with the similar length and diameter indicates 99% growth of the interfacial shear strength. Considering the possible source of errors in this level of testing, it can be claimed that there is an acceptable agreement between the experimental and the numerical results.
conclusion Multiscale damage analysis of the enhanced single fiber composite is investigated in this study. The ZnO nanowires with various geometries explored at the nanoscale are modeled as radially aligned nanowires on the carbon fiber for improving the interfacial bonding between fiber and matrix. Homogenization analysis of the proposed RVE with various geometries is implemented at the microscale to evaluate the effective material properties of the ZnO/epoxy coating layer. At mesoscale, the maximum volume fraction of the aligned ZnO on the fiber surface is explored. The Cohesive elements with bi-linear traction-separation behavior are used to simulate the interfacial bonding between the enhanced fiber and the matrix. The damage analysis of single enhanced fiber composite is performed using UMAT user subroutine at the macroscale. By comparing the results of the bare fiber with the enhanced fiber, it is concluded that the interfacial shear strength can be improved by 99% with growing ZnO on the fiber. The length and diameter effect of ZnO with common range on the interfacial strength is investigated. It is observed that smaller diameter can result in stronger interfacial bonding between fiber and matrix, which transfers the load from matrix to fiber more efficiently. Considering the effect of nanowire's length on the ZnO volume fraction in the coating layer, it is shown that shorter nanowires can result in a stronger interface. Experimental results validated the interfacial enhancement using the aligned ZnO nanowires in SFC. It is obvious that nθ = 2π, n = πd f /d ZnO . Hence, the ZnO maximum volume fraction can be re-written as Eq. A3.