A Numerical Investigation of Delamination Response of CNT/Epoxy Film Interleaved Composite

: In this study, numerical modeling through the cohesive zone theory was performed to simulate the end notch ﬂexure (ENF) test with same conditions of the experimental results of previous study that investigated the effect of a carbon nanotube (CNT)/epoxy ﬁlm in carbon ﬁber reinforced polymer (CFRP) composite through the mode II interlaminar fracture toughness of a non-interleaved, epoxy ﬁlm interleaved, CNT/epoxy ﬁlm interleaved CFRP laminate specimen. The effect of the presence of CNT/epoxy ﬁlm interleave on a composite laminate was modeled. The inﬂuence of the interleave cohesive parameters was studied to link the parameters to the material strength and energy release rate. Cohesive parameter identiﬁcation was performed by matching the initial loading and the damage evolution phase by dividing the cohesive zones into cohesive front and remaining cohesive zones. This is because, when modeling with a single cohesive zone, the critical load point that causes delamination or the curve after load drop do not match the experimental values. Results showed that the divided cohesive zone model is in good agreement with the experimental results and that there is a clear relationship between the cohesive energy of the interface and CNT/epoxy ﬁlm parameters.


Introduction
Cohesive zone modeling is based on the theory of Barenblatt and Dugdale [1,2]. The interface between the structure surfaces where crack propagation occurs is defined as the cohesive region. Fracture formation is regarded as a gradual phenomenon, wherein the separation of the surfaces involved in the crack takes place across an extended crack tip. Assuming that there is a plastic zone near the crack tip, Dugdale found the relation between plastic yield and external force. The yield strength acts across the crack. The effect of yielding is obtained as a long crack extending into the region with stress equal to the yield stress. Barenblatt allowed the stress to vary as the region in front of the crack tip deformed. The stress in the cohesive regions acts as a constraining stress holding the separating surfaces together. This stress can be referred to as cohesive stress.
There are a few studies on identifying the behavioral model based on the cohesive zone for simulating delamination of composites. Uguen et al. investigated the free-edge delamination of composite materials with the cohesive zone model. Accordingly, a comparison between different cohesive zone models was studied as well as the different parameters that influence the damage law [3]. Svensson et al. focused on the crack-tip process, which is important when studying the initiation of cracks and when crack growth has to be avoided. Various methods to measure the cohesive laws associated with the fracture processes were developed for loading in modes I and II [4]. Giuliese et al. conducted a numerical investigation of composite laminate in which Nylon 6,6 electrospun nanofibers were interleaved. The effect of nanofibers was simulated numerically by using the cohesive zone model Appl. Sci. 2022, 12, 4194 2 of 11 to investigate the difference in the fracture strengths between a nanomodified laminate interface and the corresponding non-modified material. In their study, Giuliese et al. numerically simulated double-cantilever beam (DCB) and end notch flexure (ENF) tests on CFRP composites, which has been performed in a previous study. Several configurations of nanofibrous mats have been fabricated to investigate the effect of interleaving. It was observed that the cohesive zone parameters were corresponding to the effect of nanofiber diameter, nanolayer thickness, and nanofiber orientation on the delamination behavior of the composite [5]. Hu et al. performed cohesive modeling using ABAQUS. A pre-softening zone was proposed ahead of softening zone. In this pre-softening zone, the initial stiffness and strength gradually decrease [6]. Turon et al. studied the effect of element size on force-displacement behavior in the DCB test. It shows that element size should be less than 1 mm. It is explained that this greatly affects the accuracy of the solution [7]. In addition, there are studies on the numerical modeling of delamination through cohesive zone modeling [8][9][10][11][12].
Carbon fiber reinforced polymer (CFRP) composites have high strength, stiffness, and low weight. The manufacturing techniques of CFRP composites have been developed. Thus, they substitute the existing metal materials in the aerospace, military, automotive, and sporting goods gradually.
Compared with metal materials, composite laminates are attracting attention in various industries in that they can reduce the processing cycle and enhance in-plane direction properties while making products lighter according to advances in binders and prepregs. However, in composite laminates, resin-rich areas exist between layers, resulting in inferior properties of through the thickness (or out-of-plane) direction compared to the properties of in-plane direction. Due to these characteristics, improper lamination, incomplete curing, and external impact during the process become crack seeds and delamination occurs. This failure is often considered fatal because it is often not visible from the outside. There are many studies that have aimed to improve interlaminar mechanical properties of fiber reinforced composites. Bucky paper interleaves made from carbon nanofibers improved interlaminar shear strength and mode-II interlaminar fracture toughness of composites [13]. Multi-walled carbon nanotubes (MWCNT) were grown on carbon fibers and single-walled carbon nanotubes were dispersed into matrix to enhance the interlaminar shear strength [14]. Improving the toughness of composite materials by introducing CNT has been studied by many researchers [15][16][17][18][19][20][21][22]. The interlaminar fracture toughness was improved by spraying MWCNT in a solvent, followed by evaporation and drying, and transferring the MWCNT layer to the prepreg with the inverted substrate [15].
In this work, experimental and numerical analyses of ENF tests are performed. Experimental tests described by Shin et al. [23] studied the effect of a carbon nanotube/epoxy film on the mode II interlaminar fracture toughness of CFRP laminates. The numerical simulation results were obtained through a commercial CAE analysis program, ABAQUS. Numerical modeling through the cohesive zone theory was performed to simulate the ENF test with same conditions. The influence of the carbon nanotubes (CNT)/epoxy film ply parameters was studied to link the parameters to the material strength and energy release rate. Cohesive parameter identification was performed by matching the initial loading and the damage evolution phase by dividing the cohesive zones into cohesive front and remaining cohesive zones. This study provides an effective tool to researchers and engineers who want to reinforce their laminate composite with interleaves by introducing the cohesive zone parameters obtained from the ENF tests with different concentrations of nano-reinforcements.

Materials and Experiments
The mode II interlaminar fracture toughness test described in [23] is briefly presented here. In the experiments, the laminates are made by woven carbon fiber prepreg (WSN3K, SK Chemical, Seongnam-si, Korea) and unidirectional carbon fiber prepreg (USN200A, SK Chemical). For the fabrication of CNT/epoxy film, bisphenol-A liquid epoxy resin (YD-128, Kukdo Chemical, Seoul, Korea) and urethane pre-polymer modified liquid epoxy (UME-330, Kukdo Chemical) were used as the matrix and diethyltoluenediamine (DETDA, Sigma-Aldrich, St. Louis, MO, USA) was used as the curing agent. COOH-modified MWCNT (NC3101, Nanocyl, Sambreville, Belgium) were used as the nano-reinforcement. For high content CNT/epoxy film, without solution mixing, which rapidly increases the viscosity of CNT and resin mixture, or chemical vapor deposition (CVD), which requires complicated processes, equipment, and high temperature, the CNT/epoxy mixture was first dispersed with a horn type ultrasonic probe for 12 h. After that, using a three-roll milling machine, the CNT/epoxy mixture of high viscosity according to the high content of CNT was further dispersed through shear force. For a good dispersion of CNT, it was fabricated by film casting machine through three-roll milling 5 times. An epoxy film without CNT was also fabricated for comparison. To examine the effect of interleaving of CNT/epoxy film to prevent CFRP delamination, mode II interlaminar fracture toughness values of fourteen different CFRP laminate configurations were measured (Table 1). The prepregs were cut to 200 × 200 mm and then laminated to 20 layers. After lamination of 10th layer, a 70-mm long PTFE film is placed on the front of the specimen, and the remaining 10 layers are laminated. This is for crack initiation of delamination. For specimens interleaved with epoxy film or CNT/epoxy film, the film is placed from the end of the PTFE film to the end of the specimen as shown in Figure 1. It is cured by hot-press at 140 • C for 1.5 h. More detailed process conditions and materials are described in [23]. Subsequently, it is cut to a width of 20 mm to complete the specimen production. The mode II interlaminar fracture toughness can be evaluated as follows [24]: where B is the specimen width, a is the crack length, P c is the critical load at the crack growth initialization, δ c is the critical deflection value at the loading point, and L is half the length and half the thickness of the specimen. Table 2 shows the average mode II interlaminar fracture toughness values for each type of specimen.

Numerical Modeling
The cohesive element has constitutive behavior according to the traction-separation law, as shown in Figure 2. The traction-separation law includes the initial load, the increase in stress until the onset of damage, and the decrease in stress as damage advance of the cohesive elements. When a certain level of stress or displacement is reached, the strength decreases until complete delamination of the specimen occurs. These parts can be explained by cohesive strength and cohesive energy parameters that set the properties of the cohesive region. Cohesive strength is the maximum stress required for separation. Cohesive energy is the work per unit area required to separate interface. Cohesive strength is a peak stress required for separation. Cohesive energy is separation work per unit area. The appropriate values of the parameters of the traction-separation law were chosen to describe the damage and failure of the structure.
In Figure 2a, Γ is the arbitrary path of the contour around the crack tip. The J-integral method, which represents a method for calculating the strain energy release rate along a contour, was employed to connect the cohesive theory with linear elastic fracture mechanics. (2) where is the strain energy density, is the traction vector, is the displacement vector, and is a length increment along the contour. Due to the relatively long crack length, becomes zero. Consequently, Equation (2) becomes

Numerical Modeling
The cohesive element has constitutive behavior according to the traction-separation law, as shown in Figure 2. The traction-separation law includes the initial load, the increase in stress until the onset of damage, and the decrease in stress as damage advance of the cohesive elements. When a certain level of stress or displacement is reached, the strength decreases until complete delamination of the specimen occurs. These parts can be explained by cohesive strength and cohesive energy parameters that set the properties of the cohesive region. Cohesive strength is the maximum stress required for separation. Cohesive energy is the work per unit area required to separate interface. Cohesive strength is a peak stress required for separation. Cohesive energy is separation work per unit area. The appropriate values of the parameters of the traction-separation law were chosen to describe the damage and failure of the structure.
In Figure 2a, Γ is the arbitrary path of the contour around the crack tip. The J-integral method, which represents a method for calculating the strain energy release rate along a contour, was employed to connect the cohesive theory with linear elastic fracture mechanics.
where W is the strain energy density, T is the traction vector, u is the displacement vector, and ds is a length increment along the contour. Due to the relatively long crack length, dy becomes zero. Consequently, Equation (2) becomes Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 12 * In this study, the bilinear traction-separation law is used for calculation. The initial loading and softening after damage initiation is linear in this behavior. As shown in Figure 2b, , and are the normal and tangential components of the traction vector , respectively. is the restraining stress which will increase and reach its maximum value, then begin to decrease, and eventually be reduced to zero with increasing separation.
, and are the normal and tangential components of the separation vector , respectively, at the point of the cohesive element interface. and are the initial normal and shear stiffness of the interface, respectively. Since strain energy is equivalent to the work done by external forces, is equivalent to the area under the curve, so it can be expressed by the following equation ( ). (4) Using a 2D plane strain finite element model of commercial simulation tool ABAQUS, these parameters were applied to the cohesive zone model to simulate delamination of the interleave interface of CFRP laminate by the ENF test.
To model the effect of the presence of a CNT/epoxy film interleave into an CFRP, the geometrical features of the interleave are linked to the global mechanical properties of the CFRP laminate. The mechanical properties of the unidirectional CFRP laminate and woven CFRP laminate used in simulations are presented in Table 3. The integrated square four node elements were used with a size of 1 mm. Square cohesive elements are modelled with a size of 0.5 mm, as four node two-dimensional elements (COH2D4). A maximum nominal stress criterion was applied to represent initiation and evolution of damage. A global view of modelled geometries is shown in Figure 3.  In this study, the bilinear traction-separation law is used for calculation. The initial loading and softening after damage initiation is linear in this behavior. As shown in Figure 2b, t N , and t T are the normal and tangential components of the traction vector t, respectively. t is the restraining stress which will increase and reach its maximum value, then begin to decrease, and eventually be reduced to zero with increasing separation. δ N , and δ T are the normal and tangential components of the separation vector δ, respectively, at the point of the cohesive element interface. K N and K t are the initial normal and shear stiffness of the interface, respectively. Since strain energy is equivalent to the work done by external forces, Γ is equivalent to the area under the curve, so it can be expressed by the following equation (Γ I I = G I Ic ).
Using a 2D plane strain finite element model of commercial simulation tool ABAQUS, these parameters were applied to the cohesive zone model to simulate delamination of the interleave interface of CFRP laminate by the ENF test.
To model the effect of the presence of a CNT/epoxy film interleave into an CFRP, the geometrical features of the interleave are linked to the global mechanical properties of the CFRP laminate. The mechanical properties of the unidirectional CFRP laminate and woven CFRP laminate used in simulations are presented in Table 3. The integrated square four node elements were used with a size of 1 mm. Square cohesive elements are modelled with a size of 0.5 mm, as four node two-dimensional elements (COH 2 D 4 ). A maximum nominal stress criterion was applied to represent initiation and evolution of damage. A global view of modelled geometries is shown in Figure 3.

Results and Discussion
Graphs comparing experimental and numerical values using the mechanical erties of the prepreg laminate and single cohesive zone are shown in Figure 4. Amo various cases, the result of W3 is shown as representative. Both the experimental a numerical load/displacement curves show a drop after the maximum force arrived ENF test. Delamination instantaneously overcomes the centerline of the specimen it crack propagates. It was observed that the simulated results are in good agre with the experimental results from origin to critical load. However, the numerical was inconsistent with the experimental results after the critical load was reached. case of "Simulation B", peak load is consistent with the experimental value but t gree of load drop is calculated to be relatively small. Therefore, the curve after load does not match. Conversely, in "Simulation A", the curve after load drop is match the peak load is low. As another approach in this study, two different laws of traction-separation been applied depending on the position of the cohesive zone. The length of the co zone is divided into (the length of front cohesive zone) and . The firs of the load-displacement curve, i.e., crack initiation phase, could be strongly affec the edge shape of the nano-reinforcement and PTFE film, as well as the position, dissipation, and crack initiating phenomena unlike the crack propagation phase. the cohesive zone is divided into two parts. To determine , the simulations w perimental results were conducted repeatedly. The cohesive zone parameters we with the following steps. The first step is to determine the parameters to mat maximum load. Second, identify those parameters to predict the crack propa

Results and Discussion
Graphs comparing experimental and numerical values using the mechanical properties of the prepreg laminate and single cohesive zone are shown in Figure 4. Among the various cases, the result of W3 is shown as representative. Both the experimental and the numerical load/displacement curves show a drop after the maximum force arrived in the ENF test. Delamination instantaneously overcomes the centerline of the specimen when it crack propagates. It was observed that the simulated results are in good agreement with the experimental results from origin to critical load. However, the numerical model was inconsistent with the experimental results after the critical load was reached. In the case of "Simulation B", peak load is consistent with the experimental value but the degree of load drop is calculated to be relatively small. Therefore, the curve after load drop does not match. Conversely, in "Simulation A", the curve after load drop is matched but the peak load is low.

Results and Discussion
Graphs comparing experimental and numerical values using the mechanical properties of the prepreg laminate and single cohesive zone are shown in Figure 4. Among the various cases, the result of W3 is shown as representative. Both the experimental and the numerical load/displacement curves show a drop after the maximum force arrived in the ENF test. Delamination instantaneously overcomes the centerline of the specimen when it crack propagates. It was observed that the simulated results are in good agreement with the experimental results from origin to critical load. However, the numerical model was inconsistent with the experimental results after the critical load was reached. In the case of "Simulation B", peak load is consistent with the experimental value but the degree of load drop is calculated to be relatively small. Therefore, the curve after load drop does not match. Conversely, in "Simulation A", the curve after load drop is matched but the peak load is low. As another approach in this study, two different laws of traction-separation have been applied depending on the position of the cohesive zone. The length of the cohesive zone is divided into (the length of front cohesive zone) and . The first peak of the load-displacement curve, i.e., crack initiation phase, could be strongly affected by the edge shape of the nano-reinforcement and PTFE film, as well as the position, plastic dissipation, and crack initiating phenomena unlike the crack propagation phase. Thus, the cohesive zone is divided into two parts. To determine , the simulations with experimental results were conducted repeatedly. The cohesive zone parameters were set with the following steps. The first step is to determine the parameters to match the maximum load. Second, identify those parameters to predict the crack propagation As another approach in this study, two different laws of traction-separation have been applied depending on the position of the cohesive zone. The length of the cohesive zone L is divided into L f (the length of front cohesive zone) and L − L f . The first peak of the load-displacement curve, i.e., crack initiation phase, could be strongly affected by the edge shape of the nano-reinforcement and PTFE film, as well as the position, plastic dissipation, and crack initiating phenomena unlike the crack propagation phase. Thus, the cohesive zone is divided into two parts. To determine L f , the simulations with experimental results were conducted repeatedly. The cohesive zone parameters were set with the following steps. The first step is to determine the parameters to match the maximum load. Second, identify those parameters to predict the crack propagation phase. Finally, the parameters obtained in the first step are input to elements for a length L f from the crack tip. The parameters with second step are assigned to the remaining cohesive elements.
The process of fitting the simulation results to experiments is similar to that discussed by Moroni et al. [25]. The results from the divided cohesive zone modeling are shown in Figures 5 and 6. It was observed that the simulated results are in good agreement with the experimental results beyond the critical load.      The associated cohesive parameters are presented in Table 4. In all cases, Γ II is less than Γ II since once a crack starts, less energy is required for propagation. In the case of 3, 5 wt%, both Γ II and Γ II are the highest as well as t max and t max . This is related to the fact that the CNT/epoxy film inhibits crack growth in both crack initiation and crack propagation phases when an appropriate amount of CNT is well dispersed in the epoxy. Conversely, if the amount of CNT is too large, agglomeration occurs, which leads to a decrease in stress transfer efficiency and crack propagation at the interface between resin and aggregated CNT. In addition, the Γ II and t max values of the UD laminates were higher than the Γ II and t max values of the WN laminates because of their higher moduli and strengths. By linking the geometrical features of CNT/epoxy film to the global mode II interlaminar fracture toughness properties of CFRP laminates, this represents a valuable tool for designers looking to enhance the interlaminar properties of laminates using interleaving. This finite element simulation analysis is powerful in that it can sufficiently predict the delamination behavior for non-carbon fiber or other interleaving types of reinforcement.

Conclusions
In this paper, the effect of CNT/epoxy film as an interleaving material on the ENF test was simulated numerically using the cohesive zone model. The initial loading, damage initiation with initial debonding, and damage evolution until complete separation of the CNT/epoxy film interleaved CFRP laminate under the ENF test were calculated. Even in single cohesive zone modeling, it is consistent with the increasing trend of ENF experimental results just before crack propagation. However, in order to match even after the crack initiation, a special procedure was developed that divides the cohesive zone into two and identifies each cohesive parameter. The higher G I IC , the higher the cohesive energy. This confirms that the CNT/epoxy film in which the optimal amount of CNT is dispersed inhibits the crack growth of CFRP. In addition, the finite element modeling method that can predict crack initiation and propagation more accurately is presented to researchers who want to reinforce the delamination properties of a composite laminate through interleaving.