#### Abstract

High-quality interface between propellant and insulation is the strict requirement difficult to quantify in solid rocket motor. In this study, the mixed mode delamination process of propellant and insulation interface is considered in double cantilever sandwich beam (DCSB) and single lap-joint (SLJ) test. The PPR cohesive zone model (CZM) and bilinear CZM in ABAQUS are introduced in this mixed fracture progress. In order to implement the PPR model in ABAQUS, user subroutine user element (UEL) is programmed for the novel CZM. Two simple pure mode I and mode II fracture problems are designed to check the accuracy of the UEL, and the result of verification is excellent. DCSB and SLJ test and their corresponding results are used again in the same inverse analysis with the two typical effective displacement-based and potential-based CZM. Base on the results, a series discussion and some conclusions are made. The debonding progress of the propellant and insulation interface in DCSB and SLJ test are mixed mode. The PPR CZM is prior in simulation than the bilinear CZM in ABAQUS because the PPR CZM is much more flexible with changeable traction-separation shape. The real normal and tangential displacement at damage initiation shows the unreasonable change in bilinear CZM in ABAQUS under mixed mode fracture. The PPR CZM and bilinear CZM in ABAQUS are all thickness-dependent model. The real initial stiffness and the critical displacement in the bilinear CZM and the real maximum traction in PPR model are dependent on the thickness of cohesive element. The different thickness dependence of the two model is caused by the implementation method.

#### 1. Introduction

The function of a solid rocket motor (SRM) is to deliver a thrust according to a predetermined program. SRMs are wildly used in the aerospace and missile industries to act as the launch vehicles, high-altitude escape equipment, and power device of tactical. Very strong interactions take place between propellant and insulation interface. Consequently, the properties of bonding have significant influence on performance, service life, and cost of SRMs [1]. Obviously, adhesion must satisfy many strict requirements difficult to quantify. Sufficient bonding within the entire life of the SRM should be put in the first place, for the failure interface will result in the termination of SRM’s life ahead of schedule. Many efforts have been addressed to find high-quality interface [2–4] and a series of predictive tools implemented for the assessment of the interfacial reliability [5–8].

Researchers pay more attention on finding suitable relationship to describing the debonding process, like fracture mechanics. However, the premise of a vanishingly small fracture process zone (FPZ) relative to the sample characteristic size is often not satisfied in the practical application [9]. Another alternative approach for the analysis of debonding is cohesive zone model (CZM) [10, 11]. Obviously, the fundamental issue for simulation of failure mechanisms is the characterization of cohesive interactions along the fracture surface. The mixed cohesive relationship can be classified by either nonpotential-based [12–15] or potential-based CZM [16–18]. Park and Paulino [19] pointed that effective displacement-based models may provide positive stiffness under softening conditions and the increase of separation results in the increase of corresponding traction, like the cubic polynomial traction-separation relationship proposed by Tvergaard [12]. For the potential-based CZM, several limitations are still possessed by them. The mode II fracture energy in potential-based models with polynomials is unbounded caused by the linear tangential traction [17, 20]. Similarly, the exponential-exponential potential-based models have issues on boundary conditions associated with complete failure, especially when the fracture energy of modes I and II are different [18].

Some computational methods are proposed based on the concept of CZM, like intrinsic [21] and extrinsic [22] cohesive models, extended FEM [23], microplane model [24], and virtual internal bond model [25]. The main advantages of the intrinsic cohesive model is that the model can be easily integrated with a standard finite element analysis software; thus, the intrinsic cohesive model is popular among those computational methods.

The achievement of the CZM parameters is the premise of predicting deformation and failure of the debonding interface. Most papers obtain the separation-traction relationship by assuming a priori curve combining the experimental tests performed on beam like double cantilever (DCB) and end notch flexure (ENF) [26]. If the typical experimental data, like the reaction force versus applied displacement, is identical to that simulated by inversed parameters, the input CZM parameters can be used to construct the real separation-traction relationship. On the other hand, some papers think the typical experimental data cannot represent the full field displacement [27–29]. In other words, the inversed CZM parameters must satisfy the full field displacement. However, getting full field displacement and using it in the inverse analysis are relatively complex, and the inversed results may not be much more accurate than using typical experimental data.

In order to investigate the debonding simulation of the propellant and insulation interface, intrinsic CZM combining the typical experimental data is drawn up by Zhou et al. [30] and Niu et al. [31] for mode I and mode II fracture, respectively. The double cantilever sandwich beam (DCSB) and single lap-joint (SLJ) replace the double cantilever beam (DCB) and end notched flexure (ENF) as the final test platform for the softness of propellant and insulation. Zhou et al. take the bilinear CZM in ABAQUS to obtain the mode I fracture parameters, and the mode II fracture parameters are set as 0 which is allowed in ABAQUS. For the mode II fracture, Niu et al. consider the effect of mode I and set the mode I parameters according to Zhou et al.’s results. However, for the mode I (DCSB test) and mode II (SLJ test) problems researched by Zhou et al. and Niu et al., respectively, the real model may be a mixed mode problem for the difference of materials (propellant and insulation) at the two sides of the interface. Take the mode I fracture for instance, the inconstant lateral constriction will introduce a part of tangential separation between propellant and insulation. Especially, Poisson’s ratio of the propellant always set as 0.4999 which is almost the maximum Poisson’s ratio; propellant will have a distinct lateral constriction. In addition to this, the normal separation will give projection, another part of additional tangential separation, on the deformed cohesive element. Based on the above discussion, the neglect of the mode II effect in DCSB test and the fixed mode I parameters in SLJ test is unreasonable. In other words, the cohesive element should be assigned a mixed traction-separation law in DSCB and SLJ test.

Focusing on the unreasonable factors, this paper introduces a potential-based cohesive relationship, so called PPR CZM in debonding progress simulation. In Section 2, the detailed expression of the bilinear CZM in ABAQUS and PPR CZM is introduced. Detailed implementation of the PPR CZM and its verification are conducted in Section 3. In order to compare the difference between the PPR CZM and mixed bilinear CZM in ABAQUS, same FEM models and experimental data are utilized in inverse analysis, and the corresponding result and discussion are introduced in Section 4. At last, a series conclusions are made according to the discussion.

#### 2. Cohesive Zone Model

##### 2.1. Bilinear CZM in ABAQUS

The typical bilinear traction-separation law introduced by Zhou et al. [30] is an internal built-in CZM in ABAQUS which is proposed by Davila et al. [26]. The model is defined as follows [32]: where and are the damage initiation displacement for pure mode I and pure mode II. and represent the normal and tangential cohesive strength respectively corresponding to the normal and tangential separation, and . The scalar variable represents the damage accumulated at the interface, which is 0 initially, and reaches 1, when the material is fully damaged. There are several damage evolution type in ABAQUS. Damage evolution based on effective displacement is defined as follows:

Unlike the nondimensional displacement jump vector defined in Geubelle and Baylor’s paper [15], the effective displacement is expressed as . Here, represents the Macauley bracket. The maximum value of the effective displacement is used rather than the current value to prevent healing of the interface. is the only state value that needs to be stored to record the accumulation. and are the effective displacements at complete failure and the initiation of damage, respectively.

There are four different damage initiation criterion used in ABAQUS, maximum nominal stress criterion, maximum nominal strain criterion, quadratic nominal stress criterion, and quadratic nominal strain criterion. Take the quadratic nominal strain criterion for instance. Damage is assumed to initiate when a quadratic interaction function involving the nominal strain ratios reaches a value of one. This criterion can be represented as where and represent the peak values of the nominal strain when the deformation is either purely normal to the interface or purely in the first shear direction. and are the nominal strains in the normal and first shear direction. The nominal strain in the second shear direction is neglected for the two-dimensional problems the paper focus on.

The critical normal and tangential displacement ( and ) corresponding to the fracture energy and the cohesive strength are not mentioned in ABAQUS. and are replaced by the effective displacements () at complete failure. Actually, the critical normal and tangential displacement ( and ) vary with the normal and tangential separation.

Figures 1(a) and 1(b) give the cohesive relation of the pure mode I and mode II in ABAQUS including the unloading/reloading cases represented by dotted line. The exact functional boundary of the cohesive zone model (compression in not included) is plotted in Figure 2, which is determined by the damage evolution in Equation (6).

**(a)**

**(b)**

##### 2.2. PPR Potential-Based CZM

PPR potential-based CZM is proposed by Park et al. [33] in conjunction with physical fracture parameters and consistent fracture boundary conditions. The model is expressed as where and represent the normal and tangential cohesive energy, respectively. , , , , , and are character parameters in the PPR potential-based CZM. Based on the nature of the potential function, the normal and the tangential traction are the gradient of the PPR CZM.

In addition, the unloading/reloading case can take example by the bilinear CZM in Figure 1. There are some boundaries for the novel model. The complete normal separation occurs when or in which represents the tangential conjugate final crack opening width. Similarly, the complete tangential separation occurs when or in which represents the normal conjugate final crack opening width. These boundaries define the region of the cohesive interaction as illustrated in Figure 3. As mentioned in their model, the normal and tangential conjugate final crack opening width ( and ) are the unique solution of the following nonlinear function, respectively.

**(a)**

**(b)**

Like Davila’s model, there are still two artificial stiffness indicator in the PPR CZM. The initial stiffness indicators are defined as the ratio of the of the damage initiation displacement to the critical displacements at complete failure, i.e., (, ). The character parameters and (initial slope) in the PPR CZM are expressed by the two artificial stiffness indicator as

The critical displacement can be easily obtained as follows:

In addition, the energy constants and are related to the fracture energy. When the mode I and II fracture energies are different, the energy constants and have the form

While the mode I and II fracture energies are the same, the energy constants and follow

#### 3. Computational Implantation of PPR Potential-Based CZM

There are two methods that can be used in the implementation, subroutine user element (UEL) and user material (UMAT) in ABAQUS. However, some differences lie in the two methods. UMAT is a relative simple method in programming because the program of UMAT is just a part in UEL. UMAT translates the real separation to the nominal strain and defines the constitutive model of the CZM based on the traction-separation relationship; similar handing method is used in the built-in CZM element in ABAQUS. For UEL, the traction of the user element is transformed as the equivalent nodal force in analysis. A comparative research [34] shows that UEL, though slightly more involved, is able to predict both crack ignition and large crack growth with sufficient accuracy. In addition, the UMAT consistently predicts crack initiation, but is unable to predict the crack growth accurately. Based on the above discussion, UEL is chosen for the implementation method.

##### 3.1. Basic Formulation

The formulation of cohesive interface element within the finite element method differs from standard plane elements in some aspects. First, the element has initially zero thickness. Second, only tractions across the upper and lower surfaces are involved in the balance of forces. Here, we use a linear interface element with initially zero thickness as shown in Figure 4. The weak form of the governing equation neglecting inertia and body forces and adding the contribution from the fracture interface () in the domain () can be written as where , , and are Cauchy strain tensor, displacement, and separation of the interface element, respectively. The external traction works on the boundary . is the cohesive traction along the fracture surface, the value of which will be dependent on the size of separation.

##### 3.2. Implementation in User Subroutine UEL

ABAQUS has an interface that allows users to implement special finite element in user subroutine UEL. Obviously, it is very convenient to implement user element rather than write a complete analysis code in other platform. The user element’s main contribution to the model during general analysis step is to provide residual quantity at the nodes. The residual quantity is defined as where is the external flux (caused by external traction ) and is the internal flux (due to the stress and cohesive traction along the fracture surface) at node , considering the special definition of the cohesive interface element. The residual quantity is only determined by the flux caused by cohesive traction along the fracture surface and expressed as

The local nodal displacements () as shown in Figure 5 is used to describe the separation of the interface element. First, the local displacements () should be transformed from the global displacements () by a rotational matrix (). where consists of a coordinate transformation matrix () and has the form in which is set as

Note that is the angle between the global coordinates and the local coordinates as illustrated in Figure 5. The exact definition of the and is in which where and are coordinates of the nodes before and after deformation. represents the global displacements corresponding to the local displacement in Figure 5.

The normal and tangential separation () of the surface element (see Figure 6) is defined based on the local displacements () as follows:

The relationship between the separation () and local displacements () is given as where the local displacement-separation relation matrix () is given as

Next, the separation () along the cohesive surface element is interpolated using the nodal separation and shape function. in which the shape function has the form with the linear shape function in the local coordinate which are expressed as

Combing Equations (17), (25), and (27), the separation () can be expressed by the global displacement as follows: where is global displacement-separation relation matrix.

Finally, the flux caused by cohesive traction along the fracture surface is given as in which is the real thickness of bonding surface.

In addition, the element Jacobian (stiffness matrix) is required in UEL. The stiffness matrix is defined as

##### 3.3. Simulation Verification

Based on the above analysis, a user element which used to simulate special behavior of interface cohesive zone model is implemented in UEL. Before usage, the user element should be verified by a series tests. Figure 7 gives two special tests including pure mode I (a) and pure mode II (b). The out-of-plane thickness of the problem is set as 1.0 mm. The parameters of the PPR CZM are set as , , , , , , , , and . In addition, an elastic material is used. The elastic modulus is 10000 MPa, and Poisson’s ratio is 0.3.

**(a)**

**(b)**

For pure mode I, the 2D user element is set beneath a quadrilateral plane strain element (). Displacement is applied on the top surface of the elastic element. Displacement versus time relation is plotted in Figure 8. The stress results (normal stress) and its theoretical value of user element based on element’s real separation versus displacement are illustrated in Figure 9. Firstly, while the model is compressed until the displacement is -0.05 mm, the stress is increased straightly -0.50 MPa. Then, the model is stretched until the displacement is 0 mm; the stress increased following the original path. While the displacement is increased from 0 mm to 1 mm, the user element goes through strengthen stage and reaches the maximum normal stress , and start to fracture. Then, the model is compressed; the displacement is reduced to 0 mm corresponding to the unloading case. The stress is decreased to 0 MPa and increased following the original path when the displacement increased to 1 mm again corresponding to the reloading case. Finally, the displacement increased to 3 mm, the stress decreases, and the user element is completely failure. As shown in Figure 9(a), the stress results (normal stress) are identical to the theoretical value of user element based on element’s real separation. It means that the UEL subroutine is effective in dealing with the normal traction-separation relationship.

**(a)**

**(b)**

For pure mode II, the 2D user element is set between two triangular plane strain element () as shown in Figure 7. Displacement is applied on the top and right surface of the elastic element. The corresponding displacement versus time relation is plotted in Figure 8. Unlike the pure mode I, the model is first stretched, and the displacement is increased from 0 mm to 0.5 mm. The user element is strengthen and weaken in this stage, and the maximum tangential stress is reached while displacement is 0.08 mm (see Figure 9). Then, the displacement is decreased straightly to -0.5 mm corresponding to the unloading case; the stress (tangential stress) as illustrated in Figure 9 decreased from 0.71 MPa to -0.71 MPa. After the unloading case, the model is reloaded from -0.5 mm to 0.5 mm; the corresponding stress increases following the original path. Finally, the displacement increased from 0.5 mm to 1 mm, the stress decreases, and the user element is completely failure. Obviously, the stress results coincide with the theoretical ones very well. The two simple tests mean that the UEL subroutine is effective in dealing with the mixed traction-separation relationship.

#### 4. Tests and Discussion

##### 4.1. DCSB Test

The DCSB specimen is designed in reference [30] according to the ASTM standard. The aluminum plates () are connected with the propellant plates () and insulation plates () as shown in Figure 10, and there is no debonding interface between them in order to simplify. The liner part () is used to combine the propellant and insulation. The out-of-plane thickness of the whole model is 10 mm. The model is stretched by the holes on the left. The elastic modulus of aluminum is 7000 MPa, and Poisson’s ratio is 0.33. In addition, the relaxation modulus of the propellant, , is expressed in Prony series where is the initial modulus, is the term number of the Prony series, represents the time, and and are material coefficients listed in Table 1 when the term number of the Prony series is 5. Poisson’s ratio of propellant is assumed as 0.499. Mooney-Rivlin strain energy potential model is used to describe the hyperelastic behavior of insulation. The detailed model is expressed as in which and represent the first and second deviatoric strain invariants, respectively, and is the elastic volume ratio. In addition, the material coefficients , , and are listed in Table 2.

As shown in Figure 11, the aluminum plates and propellant are modeled with four-node quadrilateral solid finite elements (CPE4R in ABAQUS), and the insulation is modeled with hybrid four-node quadrilateral solid finite elements (CPE4RH). In addition, the liner is modeled with four-node cohesive element with corresponding thickness in the viewport when using the bilinear CZM in ABAQUS. However, when UEL subroutine is used, the liner part is modeled by zero-thickness user-defined element by editing the node coordinates. The sizes of the elements and corresponding boundaries are set according to reference [30]. Displacement is applied with a constant rate of 1 mm/min. The load versus displacement was recorded until the displacement is 10 mm.

The experimental load-displacement curves for DCSB test in Figure 12 is an average curve of a series tests rather than a certain test. The inverse analysis method selected to obtain the parameters of PPR CZM and bilinear CZM in ABAQUS is Hooke-Jeeves algorithm which is also used in reference [30]. The same object function is selected in optimization with the expression

Here, and are the load from numerical and average experimental load-displacement curves. represents the displacement. is the total number of collection points. In DCSB test, the total number of collection points is 30.

Three numerical load-displacement curves for DCSB test in Figure 12 represent the bilinear CZM in reference [30], bilinear CZM in ABAQUS, and PPR CZM used in liner, respectively. The corresponding model parameters of cohesive zone models and object function are listed in Table 3. As shown in Figure 12, the maximum load in the load-displacement curves are almost the same, but the corresponding maximum normal stress are not in Table 3. The maximum normal stress of the PPR CZM is almost 5 times than the one in bilinear CZMin reference [30] and ABAQUS. Those three load-displacement curves seems to be approximation of the average experimental load-displacement curve. However, near the maximum load in load-displacement curve, the bilinear CZM in reference [30] seems to be away from the average experimental load-displacement curve. Compared with the bilinear CZM in ABAQUS and PPR CZM, the phenomenon may be caused by the ignorance of the tangential stress. As we can see in Table 3, the corresponding maximum tangential stress for bilinear CZM in ABAQUS and PPR CZM are 0.14 MPa and 0.21 MPa, respectively. This means the tangential stress has worked in contributing the numerical load-displacement curves. The object function in Table 3 shows the superiority of the PPR CZM directly.

Figure 13 gives the corresponding normal and tangential cohesive traction with respect to displacement separation in bilinear CZM in ABAQUS. Complete failure in a certain orientation (normal or tangential) is not restricted to the condition when normal or tangential separation reaches effective displacements at complete failure. A proper combination of normal or tangential separation may lead the normal or tangential traction decreased to 0 as shown in Figure 13. Normal or tangential traction’s sphere of influence is a square with a side length of . Unlike the bilinear CZM in ABAQUS, the sphere of influence for normal traction is more spacious in Figure 14. The maximum potential in PPR CZM exists in pure mode I, and the normal separation reaches the critical normal displacement . Obviously, the maximum effective normal separation for the tangential traction is much smaller than the critical normal displacement . This means the tangential traction will disappear quickly for the fast-growing normal separation in DCSB test.

**(a)**

**(b)**

##### 4.2. SLJ Test

Figure 15 gives the geometry of the DCB test for the propellant and insulation interface designed in reference [31]. Propellant and insulation plates () are still connected by the liner part (). Aluminum plates are used to act as a force transmitting device. Similarly, there is no debonding interface between aluminums and propellant (or insulation) in order to simplify. The whole model has the length of 140 mm, and the width of 16.2 mm. Unlike the DCSB test, the out-of-plane thickness of the whole model is 25 mm. As mentioned in reference [31], the material properties of the propellant is same to the one in DCSB test as well as the aluminum plates. But insulation plate has different material properties compared with DCSB test as listed in Table 2.

Finite element model and corresponding boundary conditions of the SLJ test are illustrated in Figure 16. Element types and the size of element are chosen according to reference [31]. The displacement applied at the left side will increase to 4.8 mm with a constant rate of 1 mm/min.

Average experimental load-displacement curve for SLJ test in Figure 17 is higher than numerical results before reaching the maximum load. In addition, the numerical results of different models have a similar stiffness in this stage. This phenomenon may be caused by the average effect of the average experimental load-displacement curve. At the weaken stage, PPR CZM and bilinear CZM in ABAQUS show their priority than the bilinear CZMin reference [31]. As mentioned before, the normal parameters in reference [31] is a fixed value from reference [30] in inverse analysis. The maximum normal traction is 0.31 in reference [31], and the simulated maximum normal traction is 0.28 by using the bilinear CZM in ABAQUS coincidently. This means the changed material properties of the insulation has not change the normal strength of the bonding interface much. Like the DCSB test, the maximum tangential stress in PPR CZM is 1.22 MPa which is almost 5 times than it in bilinear CZM in ABAQUS. Through much difference in the maximum tangential stress, the numerical curves and object function are similar for the three different cohesive zone model.

Figure 18 shows the normal and tangential cohesive traction with respect to displacement separation in bilinear CZM in ABAQUS in SLJ test. Unlike the bilinear CZM in ABAQUS, the tangential cohesive traction in PPR CZM increases rapidly with the increase of tangential separation at the strengthen stage in Figure 19. Obviously, the effective length of the tangential separation in PPR CZM on the normal cohesive traction is much shorter than in bilinear CZM in ABAQUS. The tangential conjugate final crack opening width is almost 10% of the critical tangential displacement as listed in Table 3. Symmetrically, the maximum PPR potential occurs with the maximum tangential separation in the case of pure mode II.

**(a)**

**(b)**

##### 4.3. Discussion

PPR CZM and bilinear CZM in ABAQUS are similar in some respects. The normal and tangential strength ( and ) parameters are utilized in their expressions. In addition, the initial stiffness is fixed by direct ( and ) or indirect (, ) input parameters. Then, bilinear CZM in ABAQUS selects an effective critical displacement to control the complete failure, and the PPR CZM uses the fracture energies ( and ) to act as the similar role. These three factors are the main governing factors of the traction-separation curve. Obviously, the two models all pay attention to those governing factors in expressions.

However, there are a series of difference between the two models. In PPR CZM, the normal (tangential) traction is dependent on the normal and tangential separation at the same time in dealing with the mixed mode separation directly. For example, the tangential separation dependence of the normal traction is considered by introducing a tangential separation-dependent coefficient in the expression of normal traction directly in Equation (5). For the bilinear CZM in ABAQUS, the mixed dependence seems relatively complicated using the damage variable . First, the damage variable will be affected by the maximum effective displacement which is controlled by normal and tangential separation. Second, the effective displacements at the initiation of damage will be changed by the ratio of normal and tangential separation. As illustrated in Figure 20, red dotted line represents the quadratic nominal strain criterion initiation criterion (the nominal strains are equal to the corresponding separation indeed), and the damage will happen when the separation path goes beyond the line like points A and B. Drawing a straight dark dotted line joining point A (or B) and origin, there will be a point of intersection with the red dotted line. The distance from the point of intersection to the origin is the effective displacements at the initiation of damage (or ). Obviously, the effective displacements at the initiation of damage will increase when tangential separation increases for displacement from the point on the red dotted line to the origin increases in Figure 20, but it will decrease when normal separation increases.

In Figure 21, normal traction with respect to the increase of the tangential separation and tangential traction with respect to the increase of the normal separation in bilinear CZM are plotted, and the parameters are from the SLJ test. The real displacements at the initiation of damage in mixed mode is increased with respect to the increased tangential separation. However, it decreases to 0 when increases from about to . Similar phenomenon happens on in Figure 21. In addition, the critical normal displacement is decreased with the increase of the tangential separation. The critical tangential displacement has the similar regulation. From Figure 21, the stiffness of the model is decreased when the tangential and normal separation are increased. Unlike the bilinear CZM in ABAQUS, the displacements at the initiation of damage and as well as the critical normal and tangential displacement and in PPR CZM never change with respect to the increase of the normal and tangential separation in Figure 22 (PPR CZM parameters are from the DCSB test). At the same time, the stiffness of the model is decreased with the increased tangential and normal separation.

**(a)**

**(b)**

**(a)**

**(b)**

Thickness dependence of the model seems never receive attention. In order to investigate the thickness dependence of the two models and the reason for the huge difference in stress results in DCSB and SLJ tests, geometry model in Figure 7(a) is selected as an simple example. Material parameters and PPR CZM parameters except thickness. The parameters for bilinear CZM in ABAQUS are , , , , and . Displacements are applied until the cohesive zone model is a complete failure. Figure 23(a) shows the strong thickness dependence of the real initial stiffness, real displacement at damage initiation, and critical displacement. Obviously, the maximum normal stress is fixed. The real initial stiffness changes because ABAQUS needs to let the nominal strain equal to the separation which is not used directly in dealing with the traction. The real initial stiffness is expressed as . At the same time, the displacement at damage initiation changes with the real initial stiffness. For the bilinear CZM, the traction will not change with the thickness for the thickness dependence of the traction transferred on the initial stiffness indirectly. In addition, ABAQUS takes the real thickness in consideration in dealing with the complete failure. The real critical displacement is the input critical displacement plus the real thickness. When changing the view to PPR CZM in Figure 23(b), initial stiffness, displacement at damage initiation, and critical displacement are not changed. However, the real maximum traction is increased with the thickness. In other words, the traction in PPR CZM is thickness dependent, and the input normal or tangential strength uses the unit thickness as the baseline. The real maximum traction is expressed as . The phenomenon is caused by the integration for the equivalent node force in Equation (31). Essentially speaking, the changed maximum traction is caused by the UEL subroutine itself. In other words, the maximum traction will not change if the UMAT subroutine is used in dealing with the PPR CZM, and the changed real stiffness may happen like the bilinear CZM in ABAQUS at the same time. Returning to the inversed results in Table 3, the real normal strength in DCSB test and the real tangential strength in SLJ test are almost same to the corresponding normal and tangential strength in other two bilinear CZMs.

**(a)**

**(b)**

As mentioned in Sections 4.1 and 4.2, the PPR CZM is excellent than the bilinear CZM in ABAQUS. In reference [30], the sensitivity analysis has been done on the bilinear CZM at the structure behavior level. However, regulations on structure behavior level cannot be transferred to different structures. From the bilinear CZM itself, as mentioned above, for the pure mode I, there are three control parameters which control three properties of the traction-separation curve, respectively, and the influence of the control parameter is independent. For the PPR CZM, as illustrated in Figure 24, , , and act the similar role as in bilinear CZM. Increasing the normal strength, the maximum normal stress will increase, and the critical displacement will reduce for the constant value of the cohesive energy Figure 24(b), and the corresponding displacement at the damage initiation will reduce. The increase of the cohesive energy will let the larger critical displacement and displacement at the damage initiation in Figure 24(c). In the above two conditions, the initial slope ( and ) will not change. Figure 24(d) shows that the increased will make the decreased initial slope and critical displacement. But the displacement at the damage initiation increases in this condition. Obviously, the three parameters can let the traction-separation curve has the same degree of freedom as the bilinear CZM, and the three parameters will not change the whole shape of the traction-separation curve. The most flexible parameter in PPR CZM, , can change the shape of the traction-separation curve at the weaken stage. This shape-control parameter lets the PPR CZM be flexible than the bilinear CZM in inverse analysis.

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

In this paper, we introduce the mixed mode cohesive zone model to realize the debonding process of the propellant and insulation interface in DCSB and SLJ test considering the asymmetric material on both sides of the interface and the insufficiency in previous research.

Typical potential-based PPR CZM and effective displacement-based bilinear CZM in ABAQUS are utilized in this paper. Detailed basic formulations for the interface element are prepared for the implementation of the PPR CZM in user subroutine UEL in ABAQUS. A series of loading and unloading/reloading cases are applied in simple pure mode I and II debonding process to demonstrate the accuracy of the UEL subroutine.

In order to make contrastive analysis, finite element model designed for the DCSB test and SLJ test, the corresponding load-displacement result and the inverse analysis method are used taking example by the previous research. The numerical results and a series discussion show that (1)The ignorance of the tangential traction in DCSB test and the fixed model parameters for mode I in SLJ test are improper in the previous research for the mixed mode debonding process existing in the two tests which is proved by the numerical results using PPR CZM and bilinear CZM in ABAQUS(2)The PPR CZM is superior than the bilinear CZM in ABAQUS based on the final object function. The advantage of the PPR CZM lies in that the model can adjust its shape by changing the shape parameter and , and the bilinear CZM in ABAQUS still keep the fixed shape(3)The normal and tangential displacement at damage initiation and the critical displacement in PPR CZM will not change in the nixed mode condition. However, for the bilinear CZM, the real normal and tangential displacement at damage initiation show the unreasonable change. For example, the real normal displacement at damage initiation increases when tangential separation is small and decreases when tangential separation is increased near the effective critical displacement. In addition, effective displacements at the initiation of damage will increase with the increased tangential separation and decrease with the increased normal separation when the bilinear CZM in ABAQUS uses the Quadratic NOMINAL strain criterion as the damage initiation in SLJ test. Similar unreasonable phenomenon will happened in DCSB test for the effective displacements at the initiation of damage(4)Bilinear CZM in ABAQUS and PPR CZM are all thickness-depended model. The real initial stiffness of the bilinear mode in ABAQUS is expressed as , and the real critical displacement is the input critical displacement plus the thickness of the cohesive element with the fixed maximum stress. The maximum stress in PPR CZM should be . The differences of the thickness dependence between the two model lie in that the method cohesive zone model implemented in ABAQUS is different. Obviously, the two methods should be extended and take the thickness into consideration like Sun et al. [35] does

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflict of interest.

#### Acknowledgments

This work is supported by the Natural Science Foundation of Jiangsu Province (BK20210435).