Peridynamic Modeling of Mode-I Delamination Growth in Double Cantilever Composite Beam Test : A Two-Dimensional Modeling Using Revised Energy-Based Failure Criteria

This study presents a two-dimensional ordinary state-based peridynamic (OSB PD) modeling of mode-I delamination growth in a double cantilever composite beam (DCB) test using revised energy-based failure criteria. The two-dimensional OSB PD composite model for DCB modeling is obtained by reformulating the previous OSB PD lamina model in x–z direction. The revised energy-based failure criteria are derived following the approach of establishing the relationship between critical bond breakage work and energy release rate. Loading increment convergence analysis and grid spacing influence study are conducted to investigate the reliability of the present modeling. The peridynamic (PD) modeling load–displacement curve and delamination growth process are then quantitatively compared with experimental results obtained from standard tests of composite DCB samples, which show good agreement between the modeling results and experimental results. The PD modeling delamination growth process damage contours are also illustrated. Finally, the influence of the revised energy-based failure criteria is investigated. The results show that the revised energy-based failure criteria improve the accuracy of the PD delamination modeling of DCB test significantly.


Introduction
Carbon fiber-reinforced polymer (CFRP) composite materials have been widely used in aerospace structures due to their high specific stiffness/strength, low coefficient of thermal expansion, and excellent fatigue resistance.In the design procedure of composite structures, testing and analysis are both needed due to the overall consideration of cost and reliability.Laminate delamination is one of the main failure modes of composite structures and, hence, has become a major concern of CFRP composite design and analysis for many years [1][2][3][4][5].Currently, the most widely used analysis method for capturing delamination initiation and growth of CFRP composites is based on the framework of finite element method (FEM), such as cohesive zone element (CZE) and virtual crack closure technique (VCCT).Although these techniques have been successfully applied to many delamination problems of CFRP composites [6][7][8], they usually have to preset the delamination growth path, which is difficult for many real engineering problems.Besides, as stated in [9][10][11][12], these conventional analysis methods based on classical continuum mechanics require that the displacement field of body be continuously differential in space.This requirement is in contradiction with inherent spatial discontinuity existing in delamination growth problems.As an alternative to conventional analysis method, Silling et al. [13][14][15] introduced the peridynamic (PD) theory of solid mechanics, which attempts to unite the mathematical modeling of continuous media, cracks, and particles within a single framework.Peridynamic theory replaces the partial differential equation of the classical theory of solid mechanics with integral or integral-differential equations, and "spontaneous" formation of fracture and damage of composites could be simulated without any prior knowledge of damage path [16][17][18][19].Peridynamic theory has been successfully applied to capture the delamination damage of CFRP composites.The interlayer delamination damage patterns of composite laminate under low velocity impact were given by Askari et al. [11] and Xu et al. [10].The delamination damages between each layer of composite laminates with notch or open hole were presented by Hu et al. [20][21][22][23]."Spontaneous" delamination damages were captured by peridynamic theory in these works, without any presetting of delamination path or refinement of the delamination front mesh.
Although the above works have proved the potential advantages of peridynamic theory in delamination modeling of CFRP composites, currently, few studies have focused on the quantitative delamination growth modeling of CFRP composite by peridynamics.Double cantilever composite beam (DCB) test provides a good benchmark example for validating an approach's basic ability in modeling mode-I delamination in composites.It is meaningful to conduct DCB modeling work to examine the reliability of an approach in modeling structural delamination problem, such as for CZE [24,25] and VCCT [26].Hu et al. [12] developed an energy-based approach to simulate delamination in composites using bond-based peridynamics.Delamination growth in double cantilever beam (DCB) test is simulated and the convex delamination front was well captured.
In the present study, a two-dimensional ordinary state-based peridynamic (OSB PD) modeling of mode-I delamination growth of CFRP composite DCB test is presented.Double cantilever beam (DCB) test is the basic experiment to measure mode-I delamination fracture toughness, G IC , of CFRP composites.This simple loading case provides a good benchmark example to quantitatively verify the ability of peridynamic theory in modeling delamination growth of CFRP composites.The two-dimensional OSB PD composite model used in the present study for DCB modeling is obtained by reformulating the previous OSB PD lamina model [27] in the x-z direction.The revised energy-based failure criteria are derived following the approach proposed by Silling et al. [15] of establishing the relationship between critical bond breakage work and energy release rate.Loading increment convergence analysis and grid spacing influence study are conducted.Comparisons of load-displacement curve and delamination growth process between PD modeling results and experimental results are shown.Delamination growth process damage contours are illustrated.Numerical results using revised and original energy-based failure criteria are also compared.

Governing Equation
The two-dimensional OSB PD composite model used in the present study for DCB modeling is obtained by reformulating the previous OSB PD lamina model (Madenci and Oterkus [27]) in x-z direction.The detailed derivation process is given in Appendix A. The governing equation of the present two-dimensional OSB PD composite model can be written as where ρ (k) is the density of material point x (k) , ü (k) is instantaneous acceleration of x (k) , as shown in Figure 1, b (k) is the external load density, and t (k)(j) and t (j)(k) are PD force density between x (k) and x (j) .The PD force density can be expressed as with and and where s (k)(j) is the stretch of bonds, and δ is the radius of the horizon zone.The direction cosines of the relative position vectors between the material points x (k) and x (j) in the undeformed and deformed states are defined as The PD dilatation θ (k) can be expressed as Figure 1.Peridynamic model notations.
The PD material parameters a, d characterize the effect of dilatation, and b, b F , b Z are associated with deformation of material points in arbitrary directions, fiber direction, and thickness direction, respectively.These parameters are related to material properties of CFRP composite and horizon radius.The detailed derivation procedures to get these PD material parameters are illustrated in Appendix A.2.These PD material parameters are related to composite material parameters as where Q 11 , Q 33 , Q 13 , and Q 55 are coefficients of stiffness matrix [Q] presented in Appendix A.1.

Revised Energy-Based Failure Criteria
Following the approach for deriving the relationship between critical bond breakage work and energy release rate by Silling et al. [15], revised two-dimensional energy-based failure criteria for mode-I delamination growth of CFRP composites are proposed.The influence of the revised energy-based failure criteria will be investigated in Section 4.3.
This approach assumes that the energy consumed by a growing delamination front equals the work required, per unit delamination front area, to separate two halves of a body across a plane (Figure 2).Suppose a plane A separates two halves of a two-dimensional body B into B + and B − .The delamination front area a is on the plane.Consider a mode-I delamination motion with velocity field on Figure 2. The total energy E absorbed by P in this motion is The assumed critical bond breakage work w IC in this motion is Therefore, For Equation ( 18), as stated by Silling et al. [15], bonds that do not cross plane A have no contributions to the integrand.However, in the present study, bonds that cross plane A are separated into three parts: bonds connecting P − and B + \P + , bonds connecting P + and B − \P − , and bonds connecting P − and P + , as shown in Figure 2. The bond breakage work between P − and P + is calculated only once compared with the original procedure in [15].Adding these three parts of the bond breakage work together gives the revised total absorbed energy E shown in Equation (18).
Assuming this total absorbed energy equals the critical energy release rate times the area of delamination front, Thus, the critical bond breakage work for mode-I delamination is related to critical energy release rate, For two-dimensional numerical modeling, the delamination front area a can be set as the discretized grid spacing dx, From the above derivation, revised energy-based failure criteria for mode-I delamination growth are proposed as Local damage at a material point is defined as the weighted ratio of the number of eliminated interactions to the total number of initial interactions of the material point with its family members [28].
The status variable for delamination damage, µ, is defined as The delamination damage between ply (n) and (n + 1) at a point can be quantified as out−of−plane_lower ). (24)

Numerical Implementation
Although the peridynamic governing equation is in dynamic form, it can still be used to solve quasi-static or static problems [29][30][31].Adaptive dynamic relaxation (ADR) method [29] is currently the most widely used method to solve quasi-static or static problems for PD.ADR method is particularly effective for solving highly nonlinear problems, including geometric and material nonlinearities [31].
According to the ADR method, Equation (1) at the n th iteration can be rewritten as where D is the fictitious diagonal density matrix and c is the damping coefficient which can be expressed by in which 1 K n is the diagonal "local" stiffness matrix, which is given as where F n i is the value of force vector F n at material point x, which includes both the peridynamic force state vector and external forces, and λ ii are the diagonal elements of D which should be large enough to avoid numerical divergence.
By utilizing central-difference explicit integration, displacements and velocities for the next time step can be obtained: and To start the iteration process, we assume that U 0 = 0 and .U 0 = 0, so the integration can be started by Due to the large computational amount of PD model, GPU-parallel computing is introduced.The PGI CUDA FORTRAN compiler, PGI/17.10Community Edition, is used for compiling.The GPU node at Cranfield University Delta HPC Cluster is applied for running the GPU-parallel program.The GPU block threads are fixed to 256, and the number of blocks depend on the total number of parallel processes [32].

Experimental Setup and Results
In order to validate the reliability of the present two-dimensional OSB PD composite model, a CFRP double cantilever beam (DCB) test is conducted.The CFRP composite DCB test setup is shown in Figure 3a.The experiment was conducted in accordance with ASTM standard D5528-13 [33].The DCB test specimens are bonded with two piano hinges at the initial delamination end.White paint marking is sprayed on the side of the specimen for visualizing the delamination front position.As shown in Figure 4, the length of the specimen is 167 mm, with 25 mm width and 4 mm thickness.The layup of the specimen is [0] 32 , with 0.125 mm/ply.Initial delamination a 0 = 50 mm is preset on the mid-plane of the specimen, between Layer 16 and 17.In this experiment, five DCB specimens are tested.The loading velocity is 1 mm/min.The resulting load-displacement curves of the specimens are presented in Figure 12.The delamination growth of each specimen is also recorded in Table 1.Modified beam theory (MBT) method from ASTM standard D5528-13 was used to analyze the experimental data, and the resulting average G IC of the CFRP composite is shown in Table 2.

Numerical Results and Discussion
The double cantilever composite beam test described in Section 3 is modeled using the two-dimensional OSB PD composite model in Section 2. The DCB specimen shown in Figure 4 is made of CFRP composite with material properties shown in Table 2.The horizon of the present two-dimensional OSB PD model is δ = 3 dx, dx is the grid spacing.The boundary conditions for PD DCB model are shown in Figure 5. Displacement boundary conditions are applied at the loading end.The loading increment for each iteration step is noted as dL.No-fail zones are set for preventing possible premature failure due to boundary effects [12].The modeling flowchart is shown in Figure 6.For each loading step, firstly, a displacement boundary condition dL is applied; secondly, ADR is used to get the static results; thirdly, the revised energy-based failure criteria derived in Section 2.2 are used to check bond failure; and finally, a new loading step is applied.The two-dimensional PD modeling deformed displacement field in thickness direction is shown in Figure 3b.
Firstly, loading increment convergence analysis for each grid spacing dx was conducted, as shown in Figures 7-10.It can be seen that with the decrease of loading increment dL, the PD modeling load-displacement curve converges for all grid spacing.As expected, the slope of the linear load-displacement part is not influenced by the loading increment, and the loading increment only influences the peak force value for a fixed grid spacing.Then, the influence of grid spacing was studied in Figure 11.The second least loading increment, dL, was used for each grid spacing.For example, dL = 0.008 mm was used for grid spacing dx = 0.25 mm.It can be seen that with the decrease of grid spacing, the linear loading part of the PD modeling load-displacement curve converges.However, PD modeling delamination growth initial point of load-displacement curve seems sensitive to grid spacing.Similar mesh sensitivity of PD modeling can be found in the work by Beckmann et al. [34], Henke et al.

Comparison of Load-Displacement Curve and Delamination Growth Process
The DCB load-displacement curve of the PD and experimental results are compared as shown in Figure 12.The load-displacement curve of PD delamination growth is shown in Figure 13.The results in Figure 12 show that the PD-based analytical result agrees well with experimental results.It is also noted that the modeling load-displacement curve of PD delamination growth shows a zigzag-shaped curve in Figure 13.This is due to the fact that when the delamination front grows for a grid spacing, the displacement-controlled load will drop.Before the next delamination front growth happens, the load will slightly increase with the increase of displacement.A similar phenomenon is also observed in the load-displacement curves obtained from experiments, as shown in Figure 12.However, due to the complex specimen manufacturing process, preset delamination condition, and experimental setup condition, the load-displacement curve obtained from experiments show significant discrepancy between different specimens.Quantitative delamination growth process is also compared between PD and experiment as presented in Table 3, and the delamination growth process damage contours are illustrated in Figure 14.In Table 3, the average delamination growth process of the five DCB test samples from 50 (initial delamination) to 65 mm (shown in Table 1) is compared with the modeling results.From the results, it noted that the present PD simulation results for the composite DCB delamination growth process agree well with the experimental results.The maximum difference between PD and experimental average results in displacement is −14.10%, and is 11.78% in load.In Figure 14, the delamination damage coefficient ϕ (16) (17) delamination is calculated using Equation (24).The delamination length growing from 50 (initial delamination) to 65 mm is vividly shown.It is worth noting that the PD delamination growth process shown in Figure 14 is "spontaneous".It does not need presetting of delamination path, nor non-physical stabilization parameters which are required for FEM modeling [26].Also, it does not require the refinement of delamination front grid spacing.

Comparison of Numerical Modeling Using Revised and Original Energy-Based Failure Criteria
In order to investigate the influence of the revised energy-based failure criteria proposed in Section 2.2, a comparison study is conducted by modeling the DCB delamination using the revised and original energy-based failure criteria.The revised energy-based failure criteria are illustrated in Equation (22).And the original energy-based failure criteria [15] for two-dimensional modeling are

Comparison of Numerical Modeling Using Revised and Original Energy-Based Failure Criteria
In order to investigate the influence of the revised energy-based failure criteria proposed in Section 2.2, a comparison study is conducted by modeling the DCB delamination using the revised and original energy-based failure criteria.The revised energy-based failure criteria are illustrated in Equation (22).And the original energy-based failure criteria [15] for two-dimensional modeling are The modeling load-displacement curve using revised and original energy-based failure criteria, as well as the experimental results, are shown in Figure 15.The comparison of maximum load using the revised and original energy-based failure criteria and the average experimental results are presented in Table 4.  From Figure 15, it is observed that the modeling load-displacement curve using revised and original energy-based failure criteria overlapped each other before delamination initiation.However, the modeling delamination initiation was delayed using the revised energy-based failure criteria, which agrees more closely with the experimental results.As shown in Table 4, differences in the maximum load of the modeling load-displacement curve in comparison with the average experimental results was reduced from −12.93% to −6.96% by using the revised energy-based failure criteria.The above results show that the revised energy-based failure criteria can significantly improve the accuracy of the PD delamination modeling of DCB test.

Conclusions
In this study, a two-dimensional ordinary state-based peridynamic (OSB PD) modeling of mode-I delamination growth in double cantilever composite beam (DCB) test is conducted.The two-dimensional OSB PD composite model for DCB modeling is obtained by reformulating the previous OSB PD lamina model in x-z direction.Additionally, revised energy-based failure criteria are proposed for modeling delamination initiation and growth.
In order to investigate the reliability of current modeling, loading increment convergence analysis and grid spacing influence study are conducted.It is shown that the PD modeling load-displacement curve converges with the decrease of loading increment.The linear loading part of PD modeling load-displacement curve converges with the decrease of grid spacing.However, the delamination growth initial point of PD modeling load-displacement curve seems sensitive to grid spacing, similar to the previous literature work.
The standard DCB test was performed with five composite DCB specimens to measure the mode-I delamination fracture toughness G IC .The load-displacement variation and delamination growths were also measured.The modeling load-displacement curve and delamination growth process using the present PD model are then compared with experimental results.The results show that the PD-based load-displacement curve agrees well with experimental results.The zigzag load-displacement curve for delamination growth was observed both in experimental results and the present PD modeling results.Well agreement of the PD modeling delamination growth progress with experimental results is achieved.In particular, the maximum difference between PD and experimental average results in displacement is −14.10%, and in load is 11.78%.Besides, the PD modeling delamination growth process damage contours are illustrated.

Figure 2 .
Figure 2. Computation of total energy absorbed by P for delamination front area a.

Table 1 .
Experimental delamination growth process observed in DCB test (DELL: delamination length).

Figure 5 .
Figure 5. Boundary conditions for PD DCB model.

Figure 7 .
Figure 7. Convergence analysis of loading increment dL for grid spacing dx = 1.0 mm, PD modeling load-displacement.

Figure 8 .
Figure 8. Convergence analysis of loading increment dL for grid spacing dx = 0.5 mm, PD modeling load-displacement.

Figure 9 .
Figure 9. Convergence analysis of loading increment dL for grid spacing dx = 0.25 mm, PD modeling load-displacement.

Figure 10 .
Figure 10.Convergence analysis of loading increment dL for grid spacing dx = 0.125 mm, PD modeling load-displacement.
[35], and Steward et al.[36].In the present study, in order to balance the computation reliability and efficiency, dx = 0.25 mm and dL = 0.008 mm are used in the following discussion.

Figure 14 .
Figure 14.PD modeling delamination growth process damage contours on the mid-plane of DCB specimen, between Layer 16 and 17.

Figure 14 .
Figure 14.PD modeling delamination growth process damage contours on the mid-plane of DCB specimen, between Layer 16 and 17.

Figure 15 .
Figure 15.Comparison of load-displacement curve obtained using revised and original bond energy-based failure criteria.

Table 2 .
Material properties of CFRP composite DCB specimen.

Table 3 .
Comparison between PD delamination growth with experimental average results (DELL: delamination length).

Table 4 .
Influence of the revised bond energy-based failure criteria on maximum load.