# The Numerical Simulation Method of the Slot Crack Groups’ Effect Analysis for Solid Rocket Motor Grain

Authors: Tianpeng Li; Xiaonan Li; Zhaolong Xuan; Yu Guo
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

This study is aimed at the stability and effect of the crack groups in the solid rocket motor (SRM) grain when it was launched at normal temperature. Based on the nearly incompressible viscoelastic finite element method, several cracks were preset in a critical location along with the dangerous point of the back slot. The singular crack elements at the tips of crack groups were established to calculate theJ-integral. With the position of the cracks, the J-integral of the various crack tips was, respectively, calculated to prejudge its stability and the group effect. Finally, the experimental measured critical J-integral JIC was compared with the numerical simulation result. The results showed that in the collinear crack groups, the enhancement effect of the main crack was caused by the nearest second crack, and the significant shielding effect of the main crack was occurred in the noncollinear crack groups. Moreover, the experimental result showed that the numerical method had high accuracy.


---

## Body

## 1. Introduction

The crack of the SRM grain was the main reason for the failure of the static test and rocket launching. The field rocket motors were required strong adaptability to the ambient temperature which was -40°C ~50°C. At the same time, the SRM grain was required to bear the combined loading which was internal pressure, temperature, and axial acceleration loading. In particular, the pressure sometimes was more than 100 atmospheres, during the SRM was launched [1, 2]. The SRM grain defects were easy to cause under serious loading conditions, such as grain crack or crack group [3, 4]. If the crack group propagation was unstable, the internal ballistic capability would be changed, and what is more, the catastrophic failure might be caused by deflagration or detonation. At present, the study of SRM grain crack was mostly limited to single crack. Actually, the detection of SRM showed that the multicracks or crack groups have commonly occurred. Therefore, it was very important to study the stability of crack groups. The propagation mode and direction of each crack in the crack group should be explored. Then, the effect of the crack group (the influence of multiple cracks on each other) would be obtained to judge the safety of SRM grain. Thus, by studying the control parameters of the crack propagation and the effect between the physical parameters of the crack group, we can get the safety effect of the crack group. This is a far-reaching research, which can broaden the horizons of researchers.The SRM is a disposable product, and its launching process is unrepeatable. In addition, the structure of the SRM has withstood the high temperature and high working internal pressure when it is launched. Therefore, it was very difficult for any effort to obtain the parameters of crack stability by experiment because of the complicated configuration and diverse loading. It is more available to establish a viscoelastic finite element method (FEM) for the analysis of the fracture parameters in SRM grain [5]. At present, with FEM, the single crack stability of the SRM grain had been discussed by many researchers [6–8], and the study on single crack propagation was also limited to the regular structure with grain crack [9]. The new methods to solve this problem also include boundary integral equation methods (BIEMs). In recent years, BIEMs have been proven to be computationally efficient and yield accurate results for the analysis of fracture problems in both 2D and 3D settings [10–16].Therefore, the effect of the crack group in SRM grain should be paid great attention. Because of the stress singularity at the crack tip, the displacement field and the first derivative of the general element were difficult to describe the real field at the crack tip, and the calculation accuracy ofJ-integral could not be improved by increasing the crack tip element. To ensure the convergence of FEM, the approximate displacement and its first derivative calculated by FEM should be arbitrarily close to the real field everywhere. The first derivative of the exact solution of the displacement field near the crack tip was unbounded; that is, the stress state of the crack tip had r−1/2 singular. A singular element with interpolation function was constructed to make the displacement of the element had r−1/2 behavior, so the stress was r−1/2 singular [17]. The displacement and stress field at the crack tip could be better described by the singular element, and calculation precision could be improved by using singular crack elements which were included in the J-integral loop.In order to study the effect of the SRM crack group, the control parameters of the crack propagation should be calculated, such as stress intensity factor [18, 19] and J-integral [3]. The J-integral of the crack tip was variated with the structure shape, crack location, and load condition of the motor grain. In this paper, a cycle-symmetric star grain configuration was taken as an example, and a method to obtain the effect of the crack group in SRM grain was proposed. The crack group of SRM grain was preset according to routine checks and service. During the whole life cycle of the SRM, the most severe load condition was the launch load condition. And the SRM was mostly used at normal temperature, and the effect of the crack group at high temperature and low temperature was not considered in this paper, so the effect of SRM grain crack group was discussed when it was launched at normal temperature [20–22]. First, the 3D element model of the SRM grain is established by using MSC.PATRAN, and the 3D singular crack elements surrounding the crack tip were developed to simulate the grain crack propagation, and all the singular crack elements were included in the 3D J-integral loop surface. Then, under combining normal temperature, internal pressure, and axial acceleration loading, the von Mises strain fields and the J-integral of each crack of the crack group were calculated by using MSC.MARC. Finally, several forms of the crack group were preset along with the SRM grain slot, and the J-integral of the variation of the crack with the position was obtained to evaluate its effect. Assuming the internal pressure value of SRM chamber and the crack cavity was the peak value.


## 2. Construction of Singular Element and Crack Element

According to fracture mechanics, the displacement and stress to the crack tip various with theJ-integral can be expressed as [23]
(1)ui=Eπ1−ν2J1/2rfiNθ,(2)σij=E2π1−ν2J1/2r1/2gijθ,where r is the distance from the crack tip to the research point (mm), E is the elastic modulus (MPa), ν is Poisson’s ratio, σiji,j=1,2,3 is the stress tension, uii=1,2,3 is the displacement tension, and fiθ and gijθ are polar functions.Equation (1) and Equation (2) show that the stress σij is obtained by calculating the first derivative of displacement ui at the crack tip. The first derivative has the singularity of r−1/2. To make the displacement and stress at the crack tip approach the true field as possible, a sort of singular element was established. The displacement of the singular element had r−1/2 behavior, so its stress which the first derivative of displacement had the singularity of r−1/2.Taking a 4-node planar quadrilateral element as an example, as shown in Figure1(a), there were two steps in the construction of a singular element. The first step is to collapse the edge to obtain the 2D singular element. It was formed by collapse node 1 and node 4 of the planar quadrilateral element as shown in Figure 1(b). The node 1 was the singular point. The second step is to surround the crack tip to obtain the crack discretization. It was formed by 16 singular elements which surrounded the crack tip as shown in Figure 1(c).Figure 1
Singular element and 2D crack element.
(a)
4-node quadrilateral element(b)
The singular element formed by degenerating the 4-node quadrilateral(c)
Building 2D crack elementThen, the displacement field of node 1 possesses the characteristic ofr, and its shape function could be expressed as
(3)N1ξ,η=1−ξ,N2ξ,η=ξ1−η,N3ξ,η=ξη.The number of singular elements around the crack tip was determined by the finite element mesh. The 2D singular elements were used to improve the calculation accuracy of planar crack fracture parameters. And all the singular crack elements were included in theJ-integral loop that could significantly improve the simulation precision [24]. The finite element meshing of the SRM grain and J-integral loop selecting are shown in Figure 2(a). The region of the crack tip was simulated with the singular crack element which was composed of 16 singular elements, and the rest was simulated with normal elements. After the stiffness matrix being assembled, the whole displacement, stress, and strain field of the SRM grain structure would be calculated.Figure 2
2D and 3D crack element meshing andJ-integral loop.
(a)
2D crack element meshing(b)
3D crack element meshing andJ-integral loopFigure2(a) shows that the loop Γ was formed by the route from bottom crack surface to top crack surface; then, the J-integral could be defined as
(4)J=∫Γωεn1−n⋅σ⋅∂u∂x1dτ,where ω is the elastic energy density, n is the outer normal of contour unit, and n1 is the component of xl direction.The 3D crack is simulated by 3D singular crack element and then calculated theJ-integral of the 3D cracks as shown in Figure 2(b). First, the 2D integral of the crack in plane (x1-x2) which was perpendicular to the crack foreside line L at pedal point O was calculated. Secondly, the 3D integral could be obtained by calculating the integral of the 2D integral integration point-by-point along the crack foreside L. The crack region was formed by the surfaces S1, S2, S3, S4, and S5; the crack region’s volume was V. An integral closed cylindrical loop surface which includes the crack element was established. Surfaces S1 and S4 were both end surface, S2 was the outer cylindrical surface, S3 was the inner cylindrical surface, S5 was the top crack surfaces, and S6 was the bottom crack surfaces.The volume domainV was composed of surface S1~S6. According to Gauss’ theorem, the area integral of the closed curved loop surface could be turned into volume integral [7]. The J-integral can be defined as
(5)J=−∫V∂∂XωεI−σ∂u∂Xρ¯dv−∫S1+S4+S5+S6σs∂u∂Xρ¯ds,where ω is the strain energy (J), n is the external normal, ρ is the unit vector and represents the crack propagation direction, ρ¯ is the weight function, and the module on the outer surface S2 is 0. On the inner surface S3, ρ¯=λlρ, between the inner and outer surfaces, ρ¯ changes smoothly between these two values, and σs is the surface tension (S1, S4, S5, and S6).Therefore, theJ-integral of each node Oi along the crack front line L can be expressed as follows [25, 26]:
(6)JOi=J∫Lλldl,where λl was the length of the crack foreside line from node Oi.


## 3. The Effect Analysis of the Slot Crack Group in the SRM Grain

According to the symmetry of the loading and configuration, as shown in Figure3(a), one-fifth of a cycle symmetric start of the SRM grain was considered for the analysis. The von Mises strain field contour chart of the SRM grain was obtained under normal temperature (20°C), internal pressure (12.5 MPa), and axial acceleration loading (20 g), as shown in Figure 3(b). The von Mises strain level of the middle zone which was shown in the A-A section was higher than the other. Hence, engineering practice showed that the axial crack group would be likely appeared at this area, no longer a single crack. In the study of the relationship between the size and accuracy of unit division, the convergence of the model is tested by the results of the maximum von Mises of propellant grain. Generally, the unit difference of 30% and the relative error of simulation results less than 5.0% are taken as the convergence criteria. As shown in Figure 1, the relationship between the engine unit division scale and the convergence of simulation results is shown. The simulation models are divided into 83152, 126186, 166112, 201658, and 268124 units, respectively. The maximum von Mises of grain obtained is 22.7%, 24.0%, 24.7%, 25.1%, and 25.3%, respectively. In order to obtain high simulation accuracy, 166112 units of division scale with relative error of 1.62% are selected.Figure 3
3D finite element model and the von Mises strain contour of the SRM being launched.
(a)
3D finite element model of the SRM(b)
The von Mises strain field contour(c)
The relationship between different element division scales and simulation convergenceDuring the whole life of the SRM, after long-term storage or low-temperature storage, the axial cracks and crack groups were easy to appear at the bottom of the SRM grain slot. Three representative crack groups were extracted to research. As shown in Figure4, the von Mises stress-strain concentration was mainly located at the bottom of the slot during the SRM being launched. The detection of cracks by solid rocket motor is complex, and the representative cases are extracted as follows. In the stress-strain concentration area, the grain crack group effects such as one or two collinear axial cracks, three collinear axial cracks, and noncollinear three cracks on the crack extension line of the main crack (the length was 100 mm, and the depth was 10 mm) tip (tip 1) were studied. The secondary crack length was b-a, the tip distance was a, and the parallel distance was h. The 3D crack element was established at each crack tip of the crack group to obtain the relationship between the J-integral of each crack tip and the position parameters a, b, and h. The J-integral criterion (JIC of this propellant grain was 0.7673 Nmm/mm2) was used to study the crack group effect.Figure 4
The slot crack groups of the SRM grain.
(a)
Double collinear cracks(b)
Three collinear cracks(c)
Three noncollinear cracksThe SRM was a composite structure consisting of various materials. The mechanical property parameters of the SRM were complicated. The propellant, insulation, and cladding were made up of polymer materials, and they had a strong time-temperature effect. The relaxation modulus of propellant, insulation, and cladding was determined from experimental measurements [26]. The expression was written as an exponential Prony series that could be written as follows:
(7)Et=E∞+∑n=1NEne−t/τn,where E∞ was the relaxation modulus at long times (t⟶∞), En was the elastic modulus of the nth Maxwell element, and τn was the relaxation time of the nth Maxwell element.The integral constitutive relation describing viscoelastic materials is as follows:(8)σt=DE0εt−εTt+D∫0tεt−s−εTt−sdEsdsds,where
(9)σ=σ11σ22σ33τ12τ23τ31T,ε=ε11ε22ε33γ12γ23γ31T,εT=αTT111000T,D=1−v1+v1−2v1v1−vv1−v000v1−v1v1−v000v1−vv1−v10000001−2v21−v0000001−2v21−v0000001−2v21−v.In order to solve the difficulty that the integral function can only be solved in the whole process of integration, the incremental constitutive relation is usually used in the finite element analysis, and the relaxation constitutive equation of formula (8) is discretized in the time domain. The relaxation modulus Et is expressed by Prony series as follows:
(10)Et=E∞+∑n=1NEne−t/τn.The time derivative of equation (10) is
(11)dEtdt=∑n=1NEnτne−t/τn.Lets=t−τ, from equation (8), the stress at t is
(12)σt=DE∞+∑n=1NEnεt−εTt+D∑n=1N−Enτn∫0tετ−εTτet−τ/τndτ.Change equation (12) as the sum of all stresses in the generalized Maxwell model
(13)σt=σ∞t+∑n=1Nσnt,where
(14)σ∞t=DE∞εt−εTt,σnt=DEnεt−εTt+−Enτn∫0tετ−εTτe−t−τ/τndτ.Time0,t is divided into N time periods with steps of h=tm−tm−1 (t0,t1, t1,t2,…, tm−1,tm,…, tN−1,tN), and the stress at tm−1 is
(15)σtm−1=σ∞tm−1+∑n=1Nσntm−1,where
(16)σ∞tm−1=DE∞εtm−1−εTtm−1,σntm−1=DEnεtm−1−εTtm−1+−Enτn∫0tm−1ετ‐εTτe−tm−1−τ/τndτ.Therefore, the stress attm can be calculated from the stress at tm−1, assuming that
(17)Δσtm=σtm−σtm−1,(18)Δεtm=εtm−εtm−1,(19)ΔεTtm=εTtm−εTtm−1.The following formula can be derived:(20)σntm=D−Enτn∫tm−1tmετ−εTτe−tm−τ/τndτ+e−h/τnσntm−1+DEn1−e−h/τnεtm−1−εTtm−1+DEnΔεtm−ΔεTtm.The total stress incrementΔσtm is
(21)Δσtm=DE∞Δεtm−ΔεTtm+∑n=1NΔσntm,where
(22)Δσntm=D−Enτn∫tm−1tmετ−εTτe−tm−τ/τndτ−1−e−h/τnσntm−1+DEn1−e−h/τnεtm−1−εTtm−1+DEnΔεtm−ΔεTtm.Suppose the strain increases linearly in time intervaltm−1,tm.
(23)εt−εTt=εtm−1−εTtm−1+t−tm+hhΔεtm−ΔεTtm.Equation (22) can be converted into
(24)Δσntm=Enτnh1−e−h/τnDΔεtm−ΔεTtm−1−e−h/τnσntm−1,αnh=1−e−h/τn,(25)βnh=τnh1−e−h/τn=αnhτnh.Then, the total stress can be expressed as(26)Δσtm=DE∞+∑n=1NβnhEnΔεtm−Δεttm−∑n=1Nαnhσntm−1.Through the above derivation, the incremental constitutive relation suitable for viscoelastic finite element analysis can be obtained.The metal case was made of steel; its elastic modulus and Poisson’s ratio were2.05×105 MPa and 0.30. The ambient temperature was 20°C. The ignition and pressurization time was 100 ms. The peak internal pressure was 12.5 MPa. The axial overload was 20 g, and Poisson’s ratio of propellant was 0.496.TheJ-integral of the main crack was affected by the second crack length of the collinear double crack, as shown in Figure 5. The J-integral of the main crack was affected by the third crack length of the collinear double crack, as shown in Figure 6. The J-integral of the main crack was affected by the distance between the main crack and the noncollinear double crack, as shown in Figure 7.Figure 5
TheJ-integral variation with second crack of double collinear cracks.Figure 6
TheJ-integral variation with third crack of three collinear cracks.Figure 7
TheJ-integral variation with the distance of three noncollinear cracks.The results showed that theJ-integral of tip 1 is 0.8096 Nmm/mm2 when the grain had a single crack (main crack), which exceeded the critical J-integral JIC of the grain, and the crack would propagate instability when the SRM being launched. Figure 5 shows that the J-integral of the main crack had decreased with the appearance of collinear crack. Therefore, the collinear crack had a shielding effect on the main crack, and the weakening extent was about 18.3%. Then, the J-integral of the main crack was less than the critical integral JIC. The main crack had become a safety crack. With the increase of the length of the second crack, the J-integral of the main crack and the two ends of the second crack tip decreased slightly. When the length of the second crack reached 25 mm, the J-integral of both the second crack and the main crack increased slowly. The increase of the length of the second crack had a certain strengthening effect on the main crack. When the length of the second crack was less than 40 mm, the enhancement was less than 9%, and the cracks would remain stable. Figure 6 shows that the second crack length of the collinear crack group was 15 mm; with the appearance of the three collinear crack groups, the J-integral of the main crack decreased by 26.1%. The effect was greater than the double crack group. As the third crack length changes, it has little effect on the J-integral of the main crack and the tip of the second crack. With the increase of the third crack length, the J-integral of the main crack and the second crack tip increased slowly. Within 50 mm of the third crack length, the enhancement amplitude was less than 2%, and each crack of the crack group was stable. Figure 7 shows that the distance between the noncollinear crack group and the main crack was less than 20 mm; all cracks of the crack group were safety cracks. As the distance (a) between the main crack and the noncollinear double crack decreased until the main crack extended 5 mm between the second crack and the third crack (a was -5 mm), the J-integral of the main crack showed a downward trend; the second crack and the third crack proximal tip (tip 2 and tip 4, respectively) had the same trend. When the distance a was 20 mm, 15 mm, 5 mm, and 0 mm, the J-integral of the main crack decreased by 10.1%, 23.0%, 26.9%, and 30.4%, respectively. Especially when the distance a was -5 mm, the J-integral of the main crack decreased by 64.7%, and it showed that the shielding phenomenon was particularly significant. As it approached the main crack, the J-integral of the distal crack tip (tip 3 and tip 5, respectively) was slightly upward trend; the enhancement amplitude was less than 3%.The results showed that theJ-integral of the main crack would be weakened when the crack group appeared at a certain distance from the front of the single crack in the SRM grain slot, and the main crack changed from unstable crack to stable crack. In the collinear crack group, the number of cracks in the crack group was proportional to the weakening amplitude, and the crack group was safe and stable. The results showed that the change of the crack length in the crack group could enhance the J-integral of the main crack in various degrees. In the collinear crack group, the second crack closest to the main crack had the greatest effect, and the effect would be weakened with the increase of the distance. The distance between the main crack and the crack group was proportional to the weakening amplitude. When the main crack tip extended the crack group area, the shielding effect was becoming more and more obvious.


## 4. Experimental Measurement and Numerical Simulation

The propellant belonged to viscoelastic material. There was no standard for the measurement of crack fracture toughness [27]. The J-integral measurement method of metal material in GB2038-91 was adopted.As shown in Figure8, the propellant unilateral crack specimen and the tensile experimental device were shown. The propellant was cut into a rectangle with a thickness of δ mm, a width of d mm, and a length of h mm, and the sizes are shown in Table 1.Figure 8
The unilateral crack specimen and tensile test device.Table 1
The sizes of propellant unilateral crack specimen.
h (mm)d (mm)δ (mm)a (mm)120.010.05.02.0, 3.0, 4.0, 5.0, 6.0Note: the test temperature was 20°C, the humidity was 69%, and the atmospheric pressure was 103.0 kPa.A unilateral crack with a length ofa was made using a blade on the side of the specimen. The J-integral of the propellant was obtained by using the method in reference [28, 29]. The formula for solving the J-integral was
(27)J=εWδd−au.In the formula,ε was the geometric shape influence factor, and W was the force-displacement curve integral. ε was calibrated by the multisample method [30–32]. According to the theory of fracture mechanics, for the J-integral under the fixed boundary displacement u was equal to
(28)J=dWdSu.S was the initial fracture ligament area of the crack body, S=ad. The formula for the ε factor obtained by the contrast formula (27) and formula (28) was
(29)ε=d−aWdWda.The multiple sample method was used to calibrate theε. The critical J-integral JIC of the propellant at different tensile rates was determined by the single sample method and formula (27). The JIC was related to the tensile rate and temperature [33, 34]. In a normal temperature environment, the JIC of the propellant at the tensile speed of 20 mm/min was measured. Five unilateral crack experimental samples with different crack lengths were prepared. The force-displacement curve was obtained by the uniaxial tensile experiment. The input energy W was obtained by integrating, and it was expressed as
(30)W=∫0ΔZPZdz.The new cast propellant specimens and the long-term stored propellant specimens were tested. TheW-a curves of 5 specimens with different crack lengths were obtained, respectively, and shown in Figure 9. The sizes of propellant unilateral crack specimen are shown in Table 2. The dW/da was obtained by fitting, and ε value was obtained by formula (29). The JIC was calculated by formula (27). The JIC of the long-term stored propellant was obtained by the experimental test to be 1.007 Nmm/mm2, and the JIC of the new cast propellant was 0.7673 Nmm/mm2.Figure 9
TheW-a curves of the new cast and the long-term stored propellant specimens.Table 2
The sizes of propellant unilateral crack specimen.
SpecimendW/da (N)ε (a=2 mm)JIC (Nmm/mm2) (average)New cast5.03380.016361.007Long-term stored4.05250.017290.7673Take the new cast propellant specimen as an example which is shown in Figure10, the 3D finite element calculation model of the propellant crack specimen was constructed. The 3D singular element was constructed at the crack tip. The finite element model at a depth of 3.0 mm was shown in the figure. The first point of the left crack tip in the figure was node 1. A total of 6 nodes form the crack front line. The radius of the J-integral was set along the crack ahead line. The critical J-integral value of the propellant crack specimen under tension at normal temperature was calculated. As shown in Figure 11, the J-integral of the crack front corresponded to the value of each node.Figure 10
The unilateral crack specimen model and 3D crack element meshing (normal temperature).Figure 11
TheJ-integral of the crack initiation (normal temperature).The crack initiation moment of the propellant crack specimen during loading was determined according to the image records in the experimental process. The critical tensile displacement and load of the crack initiation were determined by comparing the displacement-time curve in the testing equipment. It provided a reliable basis for the numerical calculation to determine the boundary conditions [35]. The average JIC of the crack tip was calculated to be 0.979 Nmm/mm2. The relative error with the test measured value was 2.8%. The results are shown in Table 3, and the analysis of the critical JIC of the propellant crack based on the viscoelastic fracture finite element method had a high accuracy.Table 3
Mechanical properties of the new cast propellant and unilateral crack results.
E∞ (MPa)Poisson’s ratioa (mm)JIC (Nmm/mm2) (average)SimulationTestRelative error1.8050.4963.00.9791.0072.8%


## 5. Conclusion

By constructing 2D and 3D crack elements of SRM grain crack, the crack groups of single crack, collinear double crack, collinear three cracks, and noncollinear three cracks at the SRM grain slot were set. TheJ-integral of each crack tip in the crack group was simulated and calculated when the SRM is launched at normal temperature. According to the relationship between the J-integral and the spatial position and size of the crack group, the effect of the crack group was discussed. The main conclusions were as follows:
(1)
In the case of a single crack, the crack was unstable and would propagate instability(2)
For collinear double cracks, the shielding effect of a crack group would appear when the second crack at a certain distance (less than 20 mm) in front of the main crack tip. At this time, both the main crack and the second crack became safety cracks. At the same time, in collinear double cracks, with the increase of the second crack length, theJ-integral of the main crack was enhanced. When the second crack length increased to 40 mm, the enhancement effect of the main crack was less than 10% after shielding, and the crack group was safe(3)
For collinear three cracks, when the length of the third crack changed to 50 mm, the enhancement effect did not exceed 2%. The enhancement effect of the third crack on the main crack should not affect its safety(4)
For noncollinear three cracks, the influence of the second crack on theJ-integral of the main crack tip was proportional to the distance. The shielding effect was gradually significant when the main crack extended the crack group from 20 mm. When the distance was -5 mm, the shielding rate was as high as 64.7%, and the crack group was more safe and stable(5)
Through the normal temperature environment, the criticalJ-integral test measured the value of the unilateral crack sample. The measured value of the critical J-integral test of the single-sided crack sample under a normal temperature environment was compared with the calculated value based on the viscoelastic fracture finite element method. The numerical simulation had high accuracy



---

*Source: 6394749-2022-11-23.xml*