Abstract
To further investigate the fracture response in propellant grain, numerical methodology is proposed to cope with crack propagation simulation especially for the mixed mode condition. The numerical discrete scheme of the propellant linear viscoelastic constitutive model is proposed, which provides a key means for the simulation of crack propagation. In order to simulate the cohesive traction distribution on the new crack surface, the extrinsic Park-Paulino-Roesler (PPR) cohesive zone model (CZM) is introduced. To let the crack propagate along any direction determined, element splitting technique and its corresponding topological operations are proposed step by step. Then, computational simulation implementation process is explained in greater detail. Typical fracture problem, single edge-notched tension test (SENT) is solved to demonstrate the efficiency and accuracy of the proposed method. In addition, double edge-notched tension test (DENT) as well as plate tension test with a slant crack is conducted to show the special fracture characters in viscoelastic solid propellant, like time dependence. Computational results reveal that the method proposed can be utilized in further fracture investigation in solid propellant combined with the experimental findings.
1. Introductions
Solid rocket motor (SRM) is widely used as the power device of rocket, missile, and sounding rocket, as well as the boost device of spacecraft launch and aircraft takeoff for its simple structure, convenient operation, and high reliability. The propellant grain is the core component of SRM. During the period from the factory to the service of SRM, various cracks, holes, blocks, scratches, and other defects may occur inside the propellant grain. These defects will produce “super” burning surface when SRM is working, which will directly or indirectly lead to the failure of the structural integrity of the drug column. The statistical results show that the number of failures caused by the damage of the structural integrity of the propellant grain accounts for 98.40% of the total SRM ignition failure. Structural integrity damage of propellant grain is the main reason for the failure of hot test or launch of solid rocket [1–3]. Because SRM has a wide range of applications in military, aerospace, and other fields, its use safety has always been a problem that cannot be ignored. To ensure the safe use of SRM, the structural integrity of the defective propellant grain must be analyzed.
The crack propagation simulation of propellant is the premise of structural integrity analysis of SRM containing defects. Since the propellant is a typical viscoelastic material, the study of crack propagation involves the intersection of viscoelastic fracture mechanics, computational solid mechanics, experimental mechanics, and other fields, and the related theories and analytical methods are still immature [4]. The above factors bring great challenges to the simulation and mechanism study of propellant crack propagation, making it one of the core and key problems in the structural integrity analysis of SRM. At present, the analysis of the integrity of SRM is still on the level of strength check of propellant grain without defect. It is of great engineering significance to carry out the simulation and mechanism research of propellant crack propagation.
Series research put efforts on the theoretical and experimental analysis of viscoelastic fracture problems. Knauss [5] derived a new differential equation to investigate the growth of a finite crack in a viscoelastic infinite plate according to Griffith’s concept of intrinsic material surface energy. Later, Schapery [6–8] proposed a theory of crack initiation and growth in viscoelastic media, including theoretical part and applications of the theory. Langlois and Gonard [9] combined the linear cumulative damage theory and the concept of failure area and derived a new law for crack propagation in solid propellant material. Experimental results reveal a strong dependence of the viscoelastic J-integral and the crack tip opening displacement for solid composite propellant [10]. Similarly, Wang et al. [11] evaluated the fracture behavior based on the J-integral which changed significantly under different tensile rates and proposed a fracture criterion considering the burning rate of the propellant.
In numerical aspect, piecewise cohesive zone model (CZM) was established to investigate the debonding, nucleation, and crack propagation in propellant [12]. Cracks are rapidly formed in the high stress area, and the redistribution of the residual stress decreased the stress level when crack propagates according to their research. However, few downtrends in total stress-strain response of the representative volume element (RVE) model means that the cracks do not propagate across the RVE model for unknown reasons. With the triangular enriched crack tip elements, Duan et al. [4, 13] investigated the deformations of crack opening and sliding displacements in the cracked viscoelastic body and obtained the energy release rate according to the enriched degree of freedom. Extended finite element method (XFEM) and finite element method (FEM) are commonly utilized in the crack propagation simulation for solid propellant. Ӧzüpek and Iyidiker [14] established a computational methodology for crack propagation analysis in solid propellant by using XFEM and CZM analysis. Han et al. [15] inserted cohesive elements between all bulk elements to simulate the crack surface in crack propagation analysis. Recently, to avoid the defect of inserting cohesive element in the whole mesh configuration, Cui [16] focused on the mode-I fracture condition in solid propellant and applied the extrinsic CZM in crack propagation simulation with FEM. Until now, computational methodology for the mixed mode crack propagation simulation of solid propellant using FEM is rarely documented.
The main objective of this study is to establish an executable computational crack propagation analysis method in detail. Viscoelastic constitutive model and its numerical implementation are introduced first in Section 2. Then, a popular CZM with the name Park-Paulino-Roesler (PPR) is presented in Section 3. In the next section, to cope with the mixed mode fracture analysis, Section 4 shows topological operations for element splitting technique. Then, a computational implementation framework is explained in detail. Numerical simulations including single edge-notched tension test (SENT), double edge-notched tension test (DENT), and plate tension test with a slant crack are investigated to prove the applicability of the proposed simulation method. Some discussions and conclusions are made in Section 7. All the analysis performed in this paper using FEM is realized in self-developed solver “CHRMULAR.”
2. Viscoelastic Constitutive Model and Its Numerical Implementation
The isotropic linear viscoelastic constitutive equation is usually expressed as follows when a constant Poisson’s ratio is considered [17]: in which , , and represent stress tensor, deviatoric stress tensor, and spherical stress tensor, respectively. and are deviatoric strain tensor and spherical strain tensor. is loading time, and is the relaxation modulus with the form of Prony series: in which represents instantaneous relaxation modulus, is the number of terms of Prony series, and and are known model parameters of the Prony series.
According to Stieltjes convolution theorem, Equations (2) and (3) can be rewritten as where and are initial deviatoric strain tensor and spherical strain tensor. The above equations will not be applied like the elastic ones directly because of the integral form in the linear viscoelastic constitutive model. In the following part, incremental form of the constitutive model will be utilized in the numerical simulation.
2.1. Deviatoric Stress Tensor Part
Considering the deviatoric part of the constitutive model, Equation (5) has the following form at .
And the form changes at .
Combining the above two equations at the two adjacent loading times, one can acquire the incremental deviatoric stress tensor.
The incremental form seems complex in the current expression. To obtain an easier formulation, three parts are designed which are expressed below:
Here
If the time interval is small enough, and will vary linearly on this time interval, and can be calculated as in which
When initial deviatoric strain tensor is ignored, a simple incremental deviatoric stress tensor can be rewritten as
2.2. Spherical Stress Tensor Part
Similarly, we can obtain the incremental spherical stress tensor.
For convenience of numerical implementation, incremental spherical stress tensor can be transformed as where
and will have a more convenient expression when time interval is small enough and can be considered to vary linearly on this time interval.
Here
When initial spherical strain tensor is ignored, a simple spherical deviatoric stress tensor can be rewritten as
According to the above numerical derivation, material matrix for the propellant with linear viscoelastic constitutive model can be written as under plane stress condition
3. Extrinsic PPR Model
As a wildly used CZM, PPR model plays a very important role in the research for solid propellant-related work and shows a good performance in simulation [18–20]. For extrinsic PPR model, the normal and tangential traction along the fracture surface are
Here, and are normal and tangential cohesive energy. and represent the opening and sliding separations. and are shape parameters. and are final crack opening and sliding separations. In addition, energy constant and have the following relationship
The final crack opening separations (, ) are controlled by the fracture energies and cohesive strength and expressed by in which, and are normal and tangential cohesive strength.
Figure 1 shows the traction response under different cohesive strength. Unlike the traditional intrinsic CZMs, traction/separation relation shows a monotonically nonincreasing curve. There is no artificial stiffness when crack propagate is much suitable for the crack propagation simulation along unknown path [21]. The final crack opening separation will decrease with the increased cohesive strength because the area between the axis and the curve is a constant value which is defined as cohesive energy.
4. Topological Operations for Element Splitting
According to Griffith’s theory that the crack will open at a plane normal to the direction of maximum stress [22], crack will propagate when maximum principal stress reaches the cohesive strength in this paper. In a previous article, crack that propagates along the element edges has been introduced in detail [16], and this part will mainly focus on the condition that crack propagates across the element. When maximum principal stress in the crack tip satisfies the criterion as state before, topological operations for element splitting will be activated following the next steps. (A)Searing all the nodes and elements in the current mesh configuration and finding the adjacent nodes and elements of node A (crack tip). In Figure 2(a), adjacent nodes and elements of node A are marked with red and light bule, respectively. For the convenience of subsequent operations, let the adjacent node direction labeled as , , , , , and . It is worth noting that and coincide with the crack surfaces which are all free edges in the current mesh. At the same time, we note the crack propagation direction with which is calculated according to maximum principal stress on node A(B)With the help of crack propagation direction and adjacent node direction, one can locate the element (marked with light green) which will be split in front of the crack tip. To avoid polygon in the current mesh when a new crack tip node D inserted after element splitting operation, adjacent element of the former split element must be split too which is marked with yellow in Figure 2(b). After the second splitting operation, two old elements (marked with light green and yellow) will be replaced by four new smaller triangular elements. Some auxiliary directions are defined in this step for the generation of new crack surface. In all node A’s new adjacent elements, connecting A and the middle point of node A’s opposite edge, adjacent element directions are labeled as , , , , and (C)In the next step, nodes on each element should be modified according to the new mesh configuration as shown in Figure 2(c). New node A’ sharing the same coordinates with node A is added to generate the new crack surface. A’s adjacent elements whose element directions lie between and should replace node A with A’ in their constitutional relations. Thus, a new crack surface generated after the participation of node A’. Considering the historical data which should be transferred during the FEM analysis, four new smaller triangular elements will inherit all the history information from their parent (elements marked with light green and yellow in Figure 2(b)).(D)At the end, to simulate the traction force along the crack surface, cohesive element is inserted using the old crack tip node A, new crack tip node D, and node A’. New crack tip node D will be used twice to form a 4-node cohesive element in Figure 2(d). To simulate the traction along the new crack surface, extrinsic PPR model as proposed in the former section will be applied on the new inserted cohesive element. Implementations of the cohesive element as well as the extrinsic PPR model can be referred in Park and Paulino’s work [23].
(a)
(b)
(c)
(d)
5. Computational Simulation Implementation
Different from the general finite element simulation calculation, the process of crack propagation simulation calculation is more complicated, and the calculation will be more time-consuming. Figure 3 shows the complete flow of propellant crack propagation simulation for readers’ reference. First, the main program loads the mesh, materials, and boundary conditions. As required, the crack tip node data that need to be tracked are marked.
Then, under the current input conditions, the corresponding displacement, stress, and strain fields are calculated. The maximum principal stress parameters of all crack tip nodes were calculated, and then, the following operations were carried out for each crack tip node. Determine whether the maximum principal stress at the crack tip node is greater than the cohesive strength. When the maximum principal stress is greater than the cohesive strength, the topological operation of crack propagation along the element edge (reference [16]) or inside the element (Section 5) is carried out. Then, the displacement, stress, and strain fields under the current mesh are recalculated. Then, the residual force of the current mesh configuration is calculated, and the convergence criterion is used to judge whether it converges. On the other hand, when the maximum principal stress is less than the cohesive strength, the convergence judgment is directly carried out.
Next, when the convergence criterion is not satisfied, the computation time step needs to be reduced, and the computation data is restored to the previous computation step. The computation was terminated when the computation time step could not be reduced any further. Conversely, when the computation converges, the current computation data are saved until the computation time reaches the total computation time.
6. Numerical Examples
6.1. Single Edge-Notched Tension Test
To verify the numerical implementation of the constitutive model, cohesive element, and element splitting technique proposed, a common single edge-notched tension test under mode-I fracture condition is utilized. In Figure 4, specimen made by solid propellant with mm and mm was designed and tested in a previous research [15]. Along its edge, a predefined crack with the length is located on the horizontal midline of the specimen. Stretch tests with displacement applied along the upper and bottom surface are conducted to investigate the fracture response of the solid propellant. Parameters in the Prony series for relaxation modulus of solid propellant are listed in Table 1 as well as adopted in the simulation. In addition, Poisson’s ratio is set as 0.49 as utilized in the reference paper [15]. In this analysis, cohesive energy and cohesive strength are 1.0 N/mm and 0.5476 MPa, respectively, in both normal and tangential direction according to reference research. Shape parameters are all set as 3.0 after parameter analysis.
Three typical mesh configurations, structured mesh (Q4), structured mesh (T3), and unstructured mesh (T3) as shown in Figure 5 are selected to check the mesh type dependency of the crack propagation simulation technique. Mesh convergence test is conducted ahead of time to avoid interference factor. Total number of elements are 4080, 15060, and 12350, respectively. 4353, 7703, and 6348 nodes are used, respectively, in those three types of mesh configuration.
Figure 6 shows the load response during the stretch test with the loading rate 10 mm/min when crack length is 15 mm. Compared with the experimental results, results shows that three different mesh types can describe the fracture response perfectly with or without element splitting technique. It is worth to mention that more elements along the crack path will provide more accurate traction on the crack surface; thus, free edges in Figures 7(b) and 7(c) seem more slippery than that in Figure 7(a). Although linear triangular element is involved, crack propagates along the right direction in this typical mode-I fracture condition which can be revealed in Figure 7(d).
(a)
(b)
(c)
(d)
Other experiments are conducted with different crack length. In Figure 8, load response is plotted when is 0.644, 0.598, 0.495, 0.441, and 0.328, respectively. Computational results are all extracted from unstructured mesh (T3) with element splitting technique. Although computational results when are not good enough, other results match the experimental results well.
With contrastive analysis of single edge-notched tension test, one can prove that numerical implementation of constitutive model, cohesive element, and element splitting technique proposed can deal with the crack simulation analysis for solid propellant effectively.
6.2. Double Edge-Notched Tension Test
A double edge-notched tension test is also employed to investigate the mixed mode fracture response under the proposed computational framework. As in Figure 9, geometry, boundary conditions, and finite element discretization of test specimen are illustrated. Refinement meshes are designed surrounding the crack tip and potential crack paths. With total 4434 nodes and 8614 triangular elements, DENT test specimen is discretized with coarse mesh and fine mesh placed in the appropriate part.
(a)
(b)
Same material parameters are adopted in this simulation. In Figure 10, load response under different loading rates is plotted. Obviously, higher loading rate will introduce higher load response directly. For viscoelastic materials, when loading time increase, relaxation modulus will decrease at the same time; thus, higher loading rate will delay the decrease of relaxation modulus compared with lower loading rate and introduce higher load response. Figure 11 shows the crack patterns of plate tension specimen when displacement is 1.5 mm. Crack propagates much deeper when mm/min as lower loading rate let the specimen show a softer texture in viscoelastic domain.
Further, under different cohesive energy and cohesive strength, parameter analysis is conducted. In Figure 12, higher cohesive energy and cohesive strength will both increase the load response. However, for changed cohesive energy, load response just various after the peak load and stay is almost unchanged before the peak load. After peak load, load response starts to show difference under different cohesive strength. Generally, higher cohesive energy consumes more fracture energy of the same specimen. On the other hand, peak load increases dramatically under higher cohesive strength. In other words, different cohesive strength will show different fracture response obviously than cohesive energy.
(a)
(b)
6.3. Plate Tension Test with a Slant Crack
Relaxation response is a typical feature for viscoelastic solid propellant. For relaxation test, a step displacement Δ is applied and keeps unchanged. Figure 13(a) shows a plate tension test specimen with width 48 mm and height 115 mm. Slant crack with angle is located on the middle of the specimen. Crack length is a fixed constant 10 mm in this simulation. Figure 13(b) shows the mesh configuration when is 45° with 4107 nodes and 8044 elements. In this section, material parameters are the same as used in the former sections.
(a)
(b)
Three types of meshes when is 30°, 60°, and 90° are selected and analyzed when step displacement Δ is 10 mm. In Figure 14, deformed mesh configurations are plotted for four typical loading times. Slant crack will be enlarged and spread out on both sides when loading time increases in the relaxation stage. Results show that cracks propagate along the horizontal direction rather than the angle of . When increases from 30° to 90°, holes introduced by crack propagation are more of a circle because the initial crack is much closer to the horizontal direction.
Figure 15 shows the effects of step displacement Δ in relaxation stage when °. When loading time increases, load response decreases rapidly at the beginning and gradually slows down. When step displacement Δ applied, strain energy is stored totally on the specimen. During the relaxation stage, strain energy is consumed on the crack propagation and deformation. Thus, larger step displacement Δ corresponding to higher strain energy will generate higher load response.
7. Conclusions
In this paper, crack propagation simulation method in viscoelastic solid propellant is investigated. Incremental method is utilized to derive the deviatoric and spherical stress tensor in each time step considering the time related integration in viscoelastic constitutive model. To take the traction along the crack surface into account, extrinsic PPR model is adopted. Topological operations for element splitting technique are presented step by step which will help the reads to have a better understanding. Computational simulation implementation is shown to provide a more executable scheme.
Numerical simulation of SENT shows a strong adaptation of the computational methodology proposed in crack propagation simulation for viscoelastic solid propellant. DENT and plate tension test with a slant crack gave convincing evidence for mixed mode crack propagation simulation. Time dependent structural response which belongs to viscoelastic material’s typical feature is demonstrated in the last two examples.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the Natural Science Foundation of Jiangsu Province (BK20210435) and the National Natural Science Foundation of China (12202496).