# Mixed Cohesive Zone Modeling of Interface Debonding between Propellant and Insulation

Authors: Huiru Cui; Weili Ma
Journal: International Journal of Aerospace Engineering
 (2022)

---

## 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.


---

## Body

## 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 Section2, 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]:
(1)Tn=σmaxΔn0Δn,Δn<0,1−DσmaxΔn0Δn,Δn≥0,Tt=1−DτmaxΔt0Δt,where Δn0 and Δt0 are the damage initiation displacement for pure mode I and pure mode II. σmax and τmax represent the normal and tangential cohesive strength respectively corresponding to the normal and tangential separation, Δn and Δt. The scalar variable D 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:
(2)D=0,Δmmax<Δm0,ΔmcΔmmax−Δm0ΔmmaxΔmc−Δm0,Δmmax≥Δm0.Unlike the nondimensional displacement jump vectorΔ¯ defined in Geubelle and Baylor’s paper [15], the effective displacement Δm is expressed as Δm=Δn2+Δt2. Here, ⋅ represents the Macauley bracket. The maximum value of the effective displacement Δmmax is used rather than the current value to prevent healing of the interface. Δmmax is the only state value that needs to be stored to record the accumulation. Δmc and Δm0 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(3)εnεn02+εtεt02=1,where εn0 and εt0 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. εn and εt 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 (Δnc and Δtc) corresponding to the fracture energy and the cohesive strength are not mentioned in ABAQUS. Δnc and Δtc are replaced by the effective displacements (Δmc) at complete failure. Actually, the critical normal and tangential displacement (Δnc and Δtc) vary with the normal and tangential separation.Figures1(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).Figure 1
Bilinear cohesive relation of the pure mode I and mode II in ABAQUS.
(a)(b)Figure 2
Exact functional boundary of the bilinear CZM in ABAQUS (compression in not included).


### 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
(4)Ψ=Δn,Δt=minϕn,ϕt+Γn1−ΔnΔncαmα+ΔnΔncm+ϕn−ϕt×Γt1−ΔtΔtcβnβ+ΔtΔtcn+ϕt−ϕn,where ϕn=∫0ΔncTnΔn,0dΔn and ϕt=∫0ΔtcTt0,ΔtdΔt represent the normal and tangential cohesive energy, respectively. α, β, Γn, Γt, m, and n 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.
(5)TnΔn,Δt=ΓnΔncm1−ΔnΔncαmα+ΔnΔncm−1−α1−ΔnΔncα−1mα+ΔnΔncm×Γt1−ΔtΔtcβnβ+ΔtΔtcn+ϕt−ϕn,(6)TtΔn,Δt=ΓtΔtcn1−ΔtΔtcβnβ+ΔtΔtcn−1−β1−ΔtΔtcβ−1nβ+ΔtΔtcn×Γn1−ΔnΔncαmα+ΔnΔncm+ϕn−ϕtΔtΔt.In addition, the unloading/reloading case can take example by the bilinear CZM in Figure1. There are some boundaries for the novel model. The complete normal separation occurs when TnΔnc,Δt=0 or TnΔn,Δtb=0 in which Δtb represents the tangential conjugate final crack opening width. Similarly, the complete tangential separation occurs when TtΔn,Δtc=0 or TtΔnb,Δt=0 in which Δnb 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 (Δnb and Δtb) are the unique solution of the following nonlinear function, respectively.
(7)fnΔnb=Γn1−ΔnbΔncmα+ΔnbΔncm+ϕn−ϕt=0,(8)ftΔtb=Γt1−ΔtbΔtcnβ+ΔtbΔtcn+ϕt−ϕn=0.Figure 3
Exact functional boundary of the PPR CZM for (a) normal traction and (b) tangential traction (compression in not included).
(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., (λn=Δn0/Δnc, λt=Δt0/Δtc). The character parameters m and n (initial slope) in the PPR CZM are expressed by the two artificial stiffness indicator as
(9)m=αα−1λn21−αλn2,n=ββ−1λt21−βλt2.The critical displacement can be easily obtained as follows:(10)Δnc=ϕnσmaxαλn1−λnα−1αm+1αmλn+1m−1,(11)Δtc=ϕtτmaxβλt1−λtβ−1βn+1βnλt+1n−1.In addition, the energy constantsΓn and Γt are related to the fracture energy. When the mode I and II fracture energies are different, the energy constants Γn and Γt have the form
(12)Γn=−ϕnϕn−ϕtϕn−ϕαmm,Γt=−ϕtϕt−ϕnϕt−ϕnβnn.While the mode I and II fracture energies are the same, the energy constantsΓn and Γt follow
(13)Γn=−ϕnαmm,Γt=βnn.


## 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]:
(1)Tn=σmaxΔn0Δn,Δn<0,1−DσmaxΔn0Δn,Δn≥0,Tt=1−DτmaxΔt0Δt,where Δn0 and Δt0 are the damage initiation displacement for pure mode I and pure mode II. σmax and τmax represent the normal and tangential cohesive strength respectively corresponding to the normal and tangential separation, Δn and Δt. The scalar variable D 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:
(2)D=0,Δmmax<Δm0,ΔmcΔmmax−Δm0ΔmmaxΔmc−Δm0,Δmmax≥Δm0.Unlike the nondimensional displacement jump vectorΔ¯ defined in Geubelle and Baylor’s paper [15], the effective displacement Δm is expressed as Δm=Δn2+Δt2. Here, ⋅ represents the Macauley bracket. The maximum value of the effective displacement Δmmax is used rather than the current value to prevent healing of the interface. Δmmax is the only state value that needs to be stored to record the accumulation. Δmc and Δm0 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(3)εnεn02+εtεt02=1,where εn0 and εt0 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. εn and εt 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 (Δnc and Δtc) corresponding to the fracture energy and the cohesive strength are not mentioned in ABAQUS. Δnc and Δtc are replaced by the effective displacements (Δmc) at complete failure. Actually, the critical normal and tangential displacement (Δnc and Δtc) vary with the normal and tangential separation.Figures1(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).Figure 1
Bilinear cohesive relation of the pure mode I and mode II in ABAQUS.
(a)(b)Figure 2
Exact functional boundary of the bilinear CZM in ABAQUS (compression in not included).


## 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
(4)Ψ=Δn,Δt=minϕn,ϕt+Γn1−ΔnΔncαmα+ΔnΔncm+ϕn−ϕt×Γt1−ΔtΔtcβnβ+ΔtΔtcn+ϕt−ϕn,where ϕn=∫0ΔncTnΔn,0dΔn and ϕt=∫0ΔtcTt0,ΔtdΔt represent the normal and tangential cohesive energy, respectively. α, β, Γn, Γt, m, and n 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.
(5)TnΔn,Δt=ΓnΔncm1−ΔnΔncαmα+ΔnΔncm−1−α1−ΔnΔncα−1mα+ΔnΔncm×Γt1−ΔtΔtcβnβ+ΔtΔtcn+ϕt−ϕn,(6)TtΔn,Δt=ΓtΔtcn1−ΔtΔtcβnβ+ΔtΔtcn−1−β1−ΔtΔtcβ−1nβ+ΔtΔtcn×Γn1−ΔnΔncαmα+ΔnΔncm+ϕn−ϕtΔtΔt.In addition, the unloading/reloading case can take example by the bilinear CZM in Figure1. There are some boundaries for the novel model. The complete normal separation occurs when TnΔnc,Δt=0 or TnΔn,Δtb=0 in which Δtb represents the tangential conjugate final crack opening width. Similarly, the complete tangential separation occurs when TtΔn,Δtc=0 or TtΔnb,Δt=0 in which Δnb 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 (Δnb and Δtb) are the unique solution of the following nonlinear function, respectively.
(7)fnΔnb=Γn1−ΔnbΔncmα+ΔnbΔncm+ϕn−ϕt=0,(8)ftΔtb=Γt1−ΔtbΔtcnβ+ΔtbΔtcn+ϕt−ϕn=0.Figure 3
Exact functional boundary of the PPR CZM for (a) normal traction and (b) tangential traction (compression in not included).
(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., (λn=Δn0/Δnc, λt=Δt0/Δtc). The character parameters m and n (initial slope) in the PPR CZM are expressed by the two artificial stiffness indicator as
(9)m=αα−1λn21−αλn2,n=ββ−1λt21−βλt2.The critical displacement can be easily obtained as follows:(10)Δnc=ϕnσmaxαλn1−λnα−1αm+1αmλn+1m−1,(11)Δtc=ϕtτmaxβλt1−λtβ−1βn+1βnλt+1n−1.In addition, the energy constantsΓn and Γt are related to the fracture energy. When the mode I and II fracture energies are different, the energy constants Γn and Γt have the form
(12)Γn=−ϕnϕn−ϕtϕn−ϕαmm,Γt=−ϕtϕt−ϕnϕt−ϕnβnn.While the mode I and II fracture energies are the same, the energy constantsΓn and Γt follow
(13)Γn=−ϕnαmm,Γt=βnn.


## 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 Figure4. The weak form of the governing equation neglecting inertia and body forces and adding the contribution from the fracture interface (Γc) in the domain (Ω) can be written as
(14)∫Ωδε:σdV+∫ΓcδΔ⋅TcdS=∫Γeδu⋅TextdS,where ε, u, and Δ are Cauchy strain tensor, displacement, and separation of the interface element, respectively. The external traction Text works on the boundary Γe. Tc is the cohesive traction along the fracture surface, the value of which will be dependent on the size of separation.Figure 4
2D cohesive zone element in local coordinates.


### 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 quantityFn is defined as
(15)Fn=Fextn−Fintn,where Fintn is the external flux (caused by external traction Text) and Fintn is the internal flux (due to the stress and cohesive traction along the fracture surface) at node n, considering the special definition of the cohesive interface element. The residual quantity Fn is only determined by the flux caused by cohesive traction along the fracture surface Fcohn and expressed as
(16)Fn=−Fcohn.The local nodal displacements (u¯) as shown in Figure 5 is used to describe the separation of the interface element. First, the local displacements (u¯) should be transformed from the global displacements (u) by a rotational matrix (R).
(17)u¯=Ru,where R consists of a coordinate transformation matrix (Λ) and has the form
(18)R=Λ0000Λ0000Λ0000Λ,in which Λ is set as
(19)Λ=cosθsinθ−sinθcosθ.Figure 5
2D cohesive zone element and displacements in the local coordinates.Note thatθ is the angle between the global coordinates and the local coordinates as illustrated in Figure 5. The exact definition of the sinθ and cosθ is
(20)sinθ=ΔyΔx2+Δy2,cosθ=ΔxΔx2+Δy2,in which
(21)x′i=xi+ui,i=1,2,3⋯8,(22)Δx=x′3+x′52‐x′1+x′72,(23)Δy=x′6+x′42‐x′8+x′22,where xi and x′i are coordinates of the nodes before and after deformation. ui represents the global displacements corresponding to the local displacement u¯i in Figure 5.The normal and tangential separation (Δ) of the surface element (see Figure 6) is defined based on the local displacements (u¯) as follows:
(24)Δ1=u¯7−u¯1,Δ2=u¯8−u¯2,Δ3=u¯5−u¯3,Δ4=u¯6−u¯4,Figure 6
The normal and tangential separation in the local coordinates.The relationship between the separation (Δ) and local displacements (u¯) is given as
(25)Δ=Lu¯,where the local displacement-separation relation matrix (L) is given as
(26)L=−100000100−100000100−101000000−10100.Next, the separation (Δ¯ξ) along the cohesive surface element is interpolated using the nodal separation and shape function.
(27)Δ¯ξ=Nu¯,in which the shape function N has the form
(28)N=N10N200N10N2,with the linear shape function in the local coordinate which are expressed as
(29)N1=1−ξ2,N2=1+ξ2.Combing Equations (17), (25), and (27), the separation (Δ¯ξ) can be expressed by the global displacement as follows:
(30)Δ¯ξ=Bu,where B=NLR is global displacement-separation relation matrix.Finally, the flux caused by cohesive traction along the fracture surfaceFcohn is given as
(31)Fcohn=∫ΓcBTTcBds=T0∫−11BTTcBdξ,in which T0 is the real thickness of bonding surface.In addition, the element Jacobian (stiffness matrix) is required in UEL. The stiffness matrixK is defined as
(32)K=−∂Fn∂u=∂Fcohn∂u=T0∫−11BT∂Tc∂uBdξ.


### 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. Figure7 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 ϕn=0.51×103N/m, ϕt=1.21×103N/m, σmax=0.51MPa, τmax=2.10MPa, α=3.7, β=3.1, λn=0.1, λn=0.05, and T0=1mm. In addition, an elastic material is used. The elastic modulus is 10000 MPa, and Poisson’s ratio is 0.3.Figure 7
Pure mode I and pure mode II test.
(a)(b)For pure mode I, the 2D user element is set beneath a quadrilateral plane strain element (1mm×1mm). 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 σmax, 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.Figure 8
Displacement versus time for pure mode I and pure mode II test.Figure 9
Stress results and its theoretical value of user element versus displacement for pure mode I (a) and pure mode II (b).
(a)(b)For pure mode II, the 2D user element is set between two triangular plane strain element (1mm×1mm×2mm) 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 τmax 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.


## 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 Figure4. The weak form of the governing equation neglecting inertia and body forces and adding the contribution from the fracture interface (Γc) in the domain (Ω) can be written as
(14)∫Ωδε:σdV+∫ΓcδΔ⋅TcdS=∫Γeδu⋅TextdS,where ε, u, and Δ are Cauchy strain tensor, displacement, and separation of the interface element, respectively. The external traction Text works on the boundary Γe. Tc is the cohesive traction along the fracture surface, the value of which will be dependent on the size of separation.Figure 4
2D cohesive zone element in local coordinates.


## 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 quantityFn is defined as
(15)Fn=Fextn−Fintn,where Fintn is the external flux (caused by external traction Text) and Fintn is the internal flux (due to the stress and cohesive traction along the fracture surface) at node n, considering the special definition of the cohesive interface element. The residual quantity Fn is only determined by the flux caused by cohesive traction along the fracture surface Fcohn and expressed as
(16)Fn=−Fcohn.The local nodal displacements (u¯) as shown in Figure 5 is used to describe the separation of the interface element. First, the local displacements (u¯) should be transformed from the global displacements (u) by a rotational matrix (R).
(17)u¯=Ru,where R consists of a coordinate transformation matrix (Λ) and has the form
(18)R=Λ0000Λ0000Λ0000Λ,in which Λ is set as
(19)Λ=cosθsinθ−sinθcosθ.Figure 5
2D cohesive zone element and displacements in the local coordinates.Note thatθ is the angle between the global coordinates and the local coordinates as illustrated in Figure 5. The exact definition of the sinθ and cosθ is
(20)sinθ=ΔyΔx2+Δy2,cosθ=ΔxΔx2+Δy2,in which
(21)x′i=xi+ui,i=1,2,3⋯8,(22)Δx=x′3+x′52‐x′1+x′72,(23)Δy=x′6+x′42‐x′8+x′22,where xi and x′i are coordinates of the nodes before and after deformation. ui represents the global displacements corresponding to the local displacement u¯i in Figure 5.The normal and tangential separation (Δ) of the surface element (see Figure 6) is defined based on the local displacements (u¯) as follows:
(24)Δ1=u¯7−u¯1,Δ2=u¯8−u¯2,Δ3=u¯5−u¯3,Δ4=u¯6−u¯4,Figure 6
The normal and tangential separation in the local coordinates.The relationship between the separation (Δ) and local displacements (u¯) is given as
(25)Δ=Lu¯,where the local displacement-separation relation matrix (L) is given as
(26)L=−100000100−100000100−101000000−10100.Next, the separation (Δ¯ξ) along the cohesive surface element is interpolated using the nodal separation and shape function.
(27)Δ¯ξ=Nu¯,in which the shape function N has the form
(28)N=N10N200N10N2,with the linear shape function in the local coordinate which are expressed as
(29)N1=1−ξ2,N2=1+ξ2.Combing Equations (17), (25), and (27), the separation (Δ¯ξ) can be expressed by the global displacement as follows:
(30)Δ¯ξ=Bu,where B=NLR is global displacement-separation relation matrix.Finally, the flux caused by cohesive traction along the fracture surfaceFcohn is given as
(31)Fcohn=∫ΓcBTTcBds=T0∫−11BTTcBdξ,in which T0 is the real thickness of bonding surface.In addition, the element Jacobian (stiffness matrix) is required in UEL. The stiffness matrixK is defined as
(32)K=−∂Fn∂u=∂Fcohn∂u=T0∫−11BT∂Tc∂uBdξ.


## 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. Figure7 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 ϕn=0.51×103N/m, ϕt=1.21×103N/m, σmax=0.51MPa, τmax=2.10MPa, α=3.7, β=3.1, λn=0.1, λn=0.05, and T0=1mm. In addition, an elastic material is used. The elastic modulus is 10000 MPa, and Poisson’s ratio is 0.3.Figure 7
Pure mode I and pure mode II test.
(a)(b)For pure mode I, the 2D user element is set beneath a quadrilateral plane strain element (1mm×1mm). 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 σmax, 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.Figure 8
Displacement versus time for pure mode I and pure mode II test.Figure 9
Stress results and its theoretical value of user element versus displacement for pure mode I (a) and pure mode II (b).
(a)(b)For pure mode II, the 2D user element is set between two triangular plane strain element (1mm×1mm×2mm) 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 τmax 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 (150mm×5mm) are connected with the propellant plates (120mm×5mm) and insulation plates (120mm×4mm) as shown in Figure 10, and there is no debonding interface between them in order to simplify. The liner part (90mm×0.2mm) 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, Et, is expressed in Prony series
(33)Et=E0−∑n=1nEn1−e−t/τnE,where E0 is the initial modulus, n is the term number of the Prony series, t represents the time, and τnE and En 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
(34)U=C10I¯1−3+C01I¯2−3+1D1Jel−12,in which I¯1 and I¯2 represent the first and second deviatoric strain invariants, respectively, and Jel is the elastic volume ratio. In addition, the material coefficients C10, C01, and D1 are listed in Table 2.Figure 10
Geometry of the DCSB test for the propellant and insulation interface.Table 1
Parameters in the Prony series for HTPB propellant stress relaxation experiment.
n012345τnE/s−−0.5921.14812.08155.579217.588En/MPa15.0008.3701.8001.5600.7340.431Table 2
The material coefficients of the hyperelastic models of the insulation.
C10/MPaC01/MPaD1/MPaDCSB test0.079660.060780.14338SLJ test0.099660.100770.11338As shown in Figure11, 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.Figure 11
Finite element model and corresponding boundary conditions of the DCSB test.The experimental load-displacement curves for DCSB test in Figure12 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
(35)minR=1N∑i=1NPΔisim−PΔiexp2.Figure 12
Average experimental and numerical load-displacement curves for DCSB test.Here,PΔisim and PΔiexp are the load from numerical and average experimental load-displacement curves. Δi represents the displacement. N is the total number of collection points. In DCSB test, the total number of collection points N is 30.Three numerical load-displacement curves for DCSB test in Figure12 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 σmax are not in Table 3. The maximum normal stress σmax 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 τmax 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.Table 3
Parameters obtained by inversion for PPR CZM, bilinear CZM in ref [30, 31], and bilinear CZM in ABAQUS in DCSB and SLJ test.
DCSB testSLJ testPPRBilinear-ref [30]Bilinear ABAQUSPPRBilinear-ref [31]Bilinear ABAQUSϕn/N⋅m‐11.970.41—0.50——ϕt/N⋅m−10.63——2.700.63—Γn/N⋅m−1-2.79——1.53——Γt/N⋅m−12.57——-2.71——σmax/MPa1.490.310.310.140.310.28τmax/MPa0.21—0.141.220.2650.25σmax,real/MPa0.30——0.03——τmax,real/MPa0.04——0.24——α3.14——2.33——β3.08——2.20——λn0.12——0.50——λt0.26——0.01——Knorm/N⋅mm−30.100.651.411.860.651.5Ktorn/N⋅mm−30.55—1.722.64e-40.902.0Δn0/mm0.360.480.222.830.480.19Δt0/mm1.47—0.080.050.290.13Δnb/mm0.50——5.66——Δtb/mm5.66——0.40——Δnc/mm3.01——5.66——Δtc/mm5.66——4.76——Δmc/mm—2.652.44—4.754.20R/N21.032.541.7738.1257.6041.26Figure13 gives the corresponding normal and tangential cohesive traction with respect to displacement separation Δn,Δt 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 Δmc. 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 Δnc. Obviously, the maximum effective normal separation Δnb for the tangential traction is much smaller than the critical normal displacement Δnc. This means the tangential traction will disappear quickly for the fast-growing normal separation in DCSB test.Figure 13
(a) Nominal and (b) tangential cohesive traction with respect to displacement separationΔn,Δt in bilinear CZM in ABAQUS in DCSB test; units: Δn, Δt, mm; Tn, Tt, MPa.
(a)(b)Figure 14
PPR potential and its gradients in DCSB test.


### 4.2. SLJ Test

Figure15 gives the geometry of the DCB test for the propellant and insulation interface designed in reference [31]. Propellant and insulation plates (40mm×5mm) are still connected by the liner part (40mm×0.2mm). 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.Figure 15
Geometry of the DCB test for the propellant and insulation interface.Finite element model and corresponding boundary conditions of the SLJ test are illustrated in Figure16. 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.Figure 16
Finite element model and corresponding boundary conditions of the SLJ test.Average experimental load-displacement curve for SLJ test in Figure17 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 σmax is 0.31 in reference [31], and the simulated maximum normal traction σmax 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 τmax 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 17
Average experimental and numerical load-displacement curves for SLJ test.Figure18 shows the normal and tangential cohesive traction with respect to displacement separation Δn,Δt 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 Δtb is almost 10% of the critical tangential displacement Δtc as listed in Table 3. Symmetrically, the maximum PPR potential occurs with the maximum tangential separation in the case of pure mode II.Figure 18
(a) Normal and (b) tangential cohesive traction with respect to displacement separationΔn,Δt in bilinear CZM in ABAQUS in SLJ test; units: Δn, Δt, mm; Tn, Tt, MPa.
(a)(b)Figure 19
PPR potential and its gradients in SLJ test.


### 4.3. Discussion

PPR CZM and bilinear CZM in ABAQUS are similar in some respects. The normal and tangential strength (σmax and τmax) parameters are utilized in their expressions. In addition, the initial stiffness is fixed by direct (Kn and Kt) or indirect (λn, λt) input parameters. Then, bilinear CZM in ABAQUS selects an effective critical displacement Δmc to control the complete failure, and the PPR CZM uses the fracture energies (ϕn and ϕt) 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 D. First, the damage variable D will be affected by the maximum effective displacement Δmmax which is controlled by normal and tangential separation. Second, the effective displacements at the initiation of damage Δm0 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 Δm0,A (or Δm0,B). Obviously, the effective displacements at the initiation of damage Δm0 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.Figure 20
Illustration of mode mix dependence ofΔm0.In Figure21, 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 Δn0 in mixed mode is increased with respect to the increased tangential separation. However, it decreases to 0 when Δt increases from about 0.7Δmc to Δmc. Similar phenomenon happens on Δt0 in Figure 21. In addition, the critical normal displacement Δnc is decreased with the increase of the tangential separation. The critical tangential displacement Δtc 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 Δn0 and Δt0 as well as the critical normal and tangential displacement Δnc and Δtc 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.Figure 21
(a) Normal traction in bilinear CZM in ABAQUS with respect to the increase of the tangential separation; (b) tangential traction in bilinear CZM in ABAQUS with respect to the increase of the normal separation.
(a)(b)Figure 22
(a) Normal traction in PPR CZM with respect to the increase of the tangential separation; (b) tangential traction in PPR CZM with respect to the increase of the normal separation.
(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 Figure7(a) is selected as an simple example. Material parameters and PPR CZM parameters except thickness. The parameters for bilinear CZM in ABAQUS are σmax=1MPa, τmax=1MPa, Δn0=1mm, Δt0=1mm, and Δmc=2mm. 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 σmax 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 Kreal=Korigin/T0. 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 Treal=Torigin⋅T0. 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.Figure 23
Illustration of thickness dependence of the normal traction of (a) bilinear CZM in ABAQUS and (b) PPR CZM.
(a)(b)As mentioned in Sections4.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, σmax, ϕn, and λn 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 (m and n) will not change. Figure 24(d) shows that the increased λn 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.Figure 24
Normal traction in PPR CZM with respect to the increase of (a)α; (b) σmax; (c) ϕn; (d) λn.
(a)(b)(c)(d)


## 4.1. DCSB Test

The DCSB specimen is designed in reference [30] according to the ASTM standard. The aluminum plates (150mm×5mm) are connected with the propellant plates (120mm×5mm) and insulation plates (120mm×4mm) as shown in Figure 10, and there is no debonding interface between them in order to simplify. The liner part (90mm×0.2mm) 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, Et, is expressed in Prony series
(33)Et=E0−∑n=1nEn1−e−t/τnE,where E0 is the initial modulus, n is the term number of the Prony series, t represents the time, and τnE and En 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
(34)U=C10I¯1−3+C01I¯2−3+1D1Jel−12,in which I¯1 and I¯2 represent the first and second deviatoric strain invariants, respectively, and Jel is the elastic volume ratio. In addition, the material coefficients C10, C01, and D1 are listed in Table 2.Figure 10
Geometry of the DCSB test for the propellant and insulation interface.Table 1
Parameters in the Prony series for HTPB propellant stress relaxation experiment.
n012345τnE/s−−0.5921.14812.08155.579217.588En/MPa15.0008.3701.8001.5600.7340.431Table 2
The material coefficients of the hyperelastic models of the insulation.
C10/MPaC01/MPaD1/MPaDCSB test0.079660.060780.14338SLJ test0.099660.100770.11338As shown in Figure11, 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.Figure 11
Finite element model and corresponding boundary conditions of the DCSB test.The experimental load-displacement curves for DCSB test in Figure12 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
(35)minR=1N∑i=1NPΔisim−PΔiexp2.Figure 12
Average experimental and numerical load-displacement curves for DCSB test.Here,PΔisim and PΔiexp are the load from numerical and average experimental load-displacement curves. Δi represents the displacement. N is the total number of collection points. In DCSB test, the total number of collection points N is 30.Three numerical load-displacement curves for DCSB test in Figure12 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 σmax are not in Table 3. The maximum normal stress σmax 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 τmax 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.Table 3
Parameters obtained by inversion for PPR CZM, bilinear CZM in ref [30, 31], and bilinear CZM in ABAQUS in DCSB and SLJ test.
DCSB testSLJ testPPRBilinear-ref [30]Bilinear ABAQUSPPRBilinear-ref [31]Bilinear ABAQUSϕn/N⋅m‐11.970.41—0.50——ϕt/N⋅m−10.63——2.700.63—Γn/N⋅m−1-2.79——1.53——Γt/N⋅m−12.57——-2.71——σmax/MPa1.490.310.310.140.310.28τmax/MPa0.21—0.141.220.2650.25σmax,real/MPa0.30——0.03——τmax,real/MPa0.04——0.24——α3.14——2.33——β3.08——2.20——λn0.12——0.50——λt0.26——0.01——Knorm/N⋅mm−30.100.651.411.860.651.5Ktorn/N⋅mm−30.55—1.722.64e-40.902.0Δn0/mm0.360.480.222.830.480.19Δt0/mm1.47—0.080.050.290.13Δnb/mm0.50——5.66——Δtb/mm5.66——0.40——Δnc/mm3.01——5.66——Δtc/mm5.66——4.76——Δmc/mm—2.652.44—4.754.20R/N21.032.541.7738.1257.6041.26Figure13 gives the corresponding normal and tangential cohesive traction with respect to displacement separation Δn,Δt 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 Δmc. 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 Δnc. Obviously, the maximum effective normal separation Δnb for the tangential traction is much smaller than the critical normal displacement Δnc. This means the tangential traction will disappear quickly for the fast-growing normal separation in DCSB test.Figure 13
(a) Nominal and (b) tangential cohesive traction with respect to displacement separationΔn,Δt in bilinear CZM in ABAQUS in DCSB test; units: Δn, Δt, mm; Tn, Tt, MPa.
(a)(b)Figure 14
PPR potential and its gradients in DCSB test.


## 4.2. SLJ Test

Figure15 gives the geometry of the DCB test for the propellant and insulation interface designed in reference [31]. Propellant and insulation plates (40mm×5mm) are still connected by the liner part (40mm×0.2mm). 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.Figure 15
Geometry of the DCB test for the propellant and insulation interface.Finite element model and corresponding boundary conditions of the SLJ test are illustrated in Figure16. 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.Figure 16
Finite element model and corresponding boundary conditions of the SLJ test.Average experimental load-displacement curve for SLJ test in Figure17 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 σmax is 0.31 in reference [31], and the simulated maximum normal traction σmax 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 τmax 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 17
Average experimental and numerical load-displacement curves for SLJ test.Figure18 shows the normal and tangential cohesive traction with respect to displacement separation Δn,Δt 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 Δtb is almost 10% of the critical tangential displacement Δtc as listed in Table 3. Symmetrically, the maximum PPR potential occurs with the maximum tangential separation in the case of pure mode II.Figure 18
(a) Normal and (b) tangential cohesive traction with respect to displacement separationΔn,Δt in bilinear CZM in ABAQUS in SLJ test; units: Δn, Δt, mm; Tn, Tt, MPa.
(a)(b)Figure 19
PPR potential and its gradients in SLJ test.


## 4.3. Discussion

PPR CZM and bilinear CZM in ABAQUS are similar in some respects. The normal and tangential strength (σmax and τmax) parameters are utilized in their expressions. In addition, the initial stiffness is fixed by direct (Kn and Kt) or indirect (λn, λt) input parameters. Then, bilinear CZM in ABAQUS selects an effective critical displacement Δmc to control the complete failure, and the PPR CZM uses the fracture energies (ϕn and ϕt) 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 D. First, the damage variable D will be affected by the maximum effective displacement Δmmax which is controlled by normal and tangential separation. Second, the effective displacements at the initiation of damage Δm0 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 Δm0,A (or Δm0,B). Obviously, the effective displacements at the initiation of damage Δm0 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.Figure 20
Illustration of mode mix dependence ofΔm0.In Figure21, 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 Δn0 in mixed mode is increased with respect to the increased tangential separation. However, it decreases to 0 when Δt increases from about 0.7Δmc to Δmc. Similar phenomenon happens on Δt0 in Figure 21. In addition, the critical normal displacement Δnc is decreased with the increase of the tangential separation. The critical tangential displacement Δtc 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 Δn0 and Δt0 as well as the critical normal and tangential displacement Δnc and Δtc 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.Figure 21
(a) Normal traction in bilinear CZM in ABAQUS with respect to the increase of the tangential separation; (b) tangential traction in bilinear CZM in ABAQUS with respect to the increase of the normal separation.
(a)(b)Figure 22
(a) Normal traction in PPR CZM with respect to the increase of the tangential separation; (b) tangential traction in PPR CZM with respect to the increase of the normal separation.
(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 Figure7(a) is selected as an simple example. Material parameters and PPR CZM parameters except thickness. The parameters for bilinear CZM in ABAQUS are σmax=1MPa, τmax=1MPa, Δn0=1mm, Δt0=1mm, and Δmc=2mm. 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 σmax 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 Kreal=Korigin/T0. 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 Treal=Torigin⋅T0. 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.Figure 23
Illustration of thickness dependence of the normal traction of (a) bilinear CZM in ABAQUS and (b) PPR CZM.
(a)(b)As mentioned in Sections4.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, σmax, ϕn, and λn 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 (m and n) will not change. Figure 24(d) shows that the increased λn 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.Figure 24
Normal traction in PPR CZM with respect to the increase of (a)α; (b) σmax; (c) ϕn; (d) λn.
(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 asKreal=Korigin/T0, 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 Treal=Torigin⋅T0. 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



---

*Source: 8772468-2022-01-19.xml*