# Precise Design of Solid Rocket Motor Heat Insulation Layer Thickness under Nonuniform Dynamic Burning Rate

Authors: Ran Wei; Futing Bao; Yang Liu; Weihua Hui
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

With the purpose of obtaining optimal designs of the heat insulating layers in solid rocket motors, we have proposed a numerical approach to compute the ideal thickness of the heat insulating layer. The proposed method is compatible with solid rocket motors that have any shape and any manner of erosion. The nonuniform dynamic burning rate is taken into consideration to achieve higher accuracy. A high-performance code is developed that uses triangular geometry as an input to allow exchanging data from any CAD platform. An improved geometric intersection algorithm is developed to generate the required sampling points, saving 35% computation time compared to its open source equivalent. Parallel computing technology is utilized to further improve the performance. All operations of the proposed approach can be executed automatically by programs, eliminating the manual work of gathering data from CAD software in the traditional approach. Validation data shows that the proposed approach saves 3.85% of the mass compared to the ordinary design approach. Performance profiling shows that the implemented code operates within 5 seconds, which is much faster than the unoptimized open source version.


---

## Body

## 1. Introduction

The heat insulating layer (HIL) is an important component of a solid rocket motor (SRM) used to prevent the high-temperature gas from attacking the case. The burning and recession of the SRM grain and HIL during the working process of an SRM can be summarized by Figure1. It can be seen that as the grain recesses the HIL gets gradually exposed. Since the exposed part of the HIL is continuously attacked by the high-temperature gas and condensed particles, the HIL is eroded slowly during the procedure. The protective effect lasts until the residual HIL becomes too thin to resist heat. Ideally, the thickness of the residual HIL after combustion will exactly meet the minimum safety requirements. Extra heat insulating material introduces unnecessary weight, while insufficient heat insulating material brings the risk of the case being burnt out. As a result, SRM designers have to predict the erosion at each location to decide the optimal design of HILs, which are thick enough to resist heat during the process while being as light as possible to reduce the weight of SRMs.Figure 1
Burning and recession of the grain and HIL.This article is addressed to develop an approach to provide such optimal designs of HILs. In order to compute the ideal thickness, two basic data sets are required: (a) the erosion characteristics of the specific heat insulating material in SRMs and (b) the duration that each point on the HIL is exposed to the gas stream and particles. In both of the fields, there have been a number of basic studies. However, we have noticed that few studies have been aimed at integrating the work across the two fields to establish an upper-level tool. In quite a number of engineering practices, constant erosion and combustion rate presumption, manual geometric analysis, and CAD measurement are utilized to estimate the optimal HIL thickness. Such practices introduce nonnegligible error and therefore impels the designers to choose larger safety factors for HILs to compensate for the possible error. Besides, relying on CAD platforms to gather computation input usually introduces trivial manual measurement work. The absence of automatic universal upper-level designing tools drives us to develop the approach presented in this article.Previous studies focused on the erosion of HILs were comprehensively summarized by Xu et al. [1]. There are simple heat ablation approaches [2, 3], thermochemical ablation methods [4, 5], and methods which take condensed particles into consideration [6–8]. Each of these works applies to specific thermochemical ablation environments and provides an accurate prediction of the local erosion rate. In the sense of predicting the eroded thickness, despite the complex physical and mathematical models behind these erosion models, any erosion model can be characterized by a time-dependent function describing the recession rate of the HIL. In this article, an erosion model is summarized by a recession rate function which takes the material properties, the position, and the time as its arguments. Such encapsulation evades the full erosion simulation of the SRM by decoupling the HIL designing from erosion simulation, so the designers are free to integrate whatever erosion model into the presented approach.The exposure duration of a particular location depends on the transient burning rate and the distance the flame traveled from its initial position to this location. Since the burning rate of a propellant varies as the inner pressure fluctuates, and that the shortest path for the flame to propagate is not always a straight line, computing the exact exposure duration becomes a knotty work. In this article, we have proposed a numerical technology, which utilizes SRM grain burnback models to find the shortest path of the flame. An integration-based method is then developed to compute the exposure duration.Earlier studies about the burnback model can be divided into two categories: (a) the uniform models [9, 10], which assume the transient burning rate to be consistent across the entire grain, and (b) the nonuniform models [11–14], which take the spatial variation of the burning rate (caused by erosive burning or multiple types of propellants) into consideration. Generally, uniform burnback models handle most SRMs at enough accuracy, while the nonuniform ones provide better versatility at the cost of much higher computation consumption. In this article, we have developed a high-performance implementation of the uniform burnback model, while the level set method [12] is used as the nonuniform burnback. The work presented by [12] only produces burning-area data with respect to flame propagation distance, so we have enhanced the implementation to predict the exposure duration data as well.To sum up, we have developed a general approach, which bridges the previous studies on HIL erosion and grain burnback, to automatically produce optimal designs of the HILs. The time-dependent, possibly spatially nonuniformly distributed burning rate is taken into consideration to improve the accuracy compared to the common engineering practice. Both the erosion model and the burnback model are replaceable to fit different SRMs. Besides, we have also presented a high-performance implementation of the proposed approach, which makes it possible to embed the approach into higher-level optimizing frameworks. In Section2 of this article, we firstly build a universal mathematical model covering all relevant physical processes. Then, in Section 3, the details about the high-performance code implementation is provided. Finally, in Section 4, we have validated the established code and demonstrated the improvements of our approach in comparison to the traditional method.


## 2. Mathematical Model

### 2.1. Universal Model

Our study is based on the following preconditions and properties of the SRM:P.1.
The SRM is discretized to a triangular model with enough accuracy. Each of the input triangles is tagged with its associated information
P.2.
The grain is case bonded. If not, the entire HIL would be exposed to the high-temperature gas from the beginning of burning, and the simulation would be unnecessary
P.3.
There is no hole or flame-retardant material inside the grain
P.4.
The thickness of the HIL on the same radial section is uniform. Otherwise, manufacturing of the HIL and the grain would become difficult
P.5.
The shape change of the HIL caused by erosion is negligible compared to its overall scaleAccording to the associated information, we define the set of all triangles on the initial burning surface asB, and the set of the triangles on the inner surface of HIL as H.LettingTEx→ be the duration before the flame arriving point x→, we have
(1)Tx→=∫Px→1rt,s→ds,where Px→ is the shortest flame traveling path from B to x→ and rt,s→ is the transient burning rate at the burning location on Px→. Since the burning rate has spatial dependency, the shortest path Px→ may not be a straight line in a general case, where nonuniform burnback models have to be integrated. For the SRMs which can be approximately handled by the uniform burnback model, Px→ is always the shortest line segment connecting B and x→, so the solution of equation (1) can be greatly simplified (see the discussion in Section 2.2).The burning ratert,x→ is decided by the widely known exponential burning rate formula [15]:
(2)rt,x→=Ax→pt,x→nx→,where Ax→ and nx→ are decided by the local propellant material property and initial temperature and pt,x→ is the current gas pressure near x→. Usually, the values of Ax→ and nx→ can be determined by burning rate experiments and are constant among different SRMS, while pt,x→ has to be computed by steady or transient inner ballistic models.For the HIL, the ablation process starts the moment the HIL is exposed to the gas. Lettingte be the elapsed time after exposure, the erosion model at x→ can be summarized by the following time-dependent function:
(3)RRte=MEHtcdtte,Htadtte,Mcte,HM,te,where RRte denotes the recession rate at time te, Htcdt and Htadt denote the heat transfer caused by conduction and advection, respectively, Mc denotes the mechanical effect of condensed particles, and HM stands for the material property of the HIL. The last parameter te implies that function ME has a time-memory property, i.e., ME may output different values at different moments even if Htcdt, Htadt, and Mc remain identical.LetTend be the moment when the SRM stops working. The thickness of the eroded HIL can be computed as follows:
(4)Ex→=∫0Tend−TEx→RRtedte.In order to prevent the overheating of SRM cases, the thickness of the residual HIL atTend must not be too small. The ideal residual thickness can be decided according to the material properties of the SRM. Letting L0 be the minimum allowed residual thickness, we have the following to calculate the optimal thickness of HIL:
(5)Lx→=L0+Ex→.


### 2.2. Model Specialization

Section2.1 presents the universal model which theoretically applies to any SRM. Practically, the performance of most SRMs can be approximately simulated by the uniform burning rate model [9, 10], which assumes that the change in pt,x→ with the position is negligible and thus uses a consistent transient burning rate across the entire grain.By removing the spatial dependency from the burning rate formula, equation (1) can be simplified as the following form:
(6)Dx→,B=∫0Tx→Aptndt,where Dx→,B is the length of the shortest line segment connecting B to x→. There have been some high-performance minimum distance algorithms, including the Havoc algorithm [16] employed in [10], the DiFi algorithm [17] employed in [12], and the Axis-Aligned Bounding Box (AABB) tree method [18]. Both Havoc and DiFi rely on the Voronoi diagram [19] to query the minimum distance. In this article, only the points on H are of our concern, so there are fewer points to process compared to full burnback simulations [10, 12]. Our tests show that the AABB tree method outperforms the Voronoi-type methods when the number of points to be queried is relatively small.


## 2.1. Universal Model

Our study is based on the following preconditions and properties of the SRM:P.1.
The SRM is discretized to a triangular model with enough accuracy. Each of the input triangles is tagged with its associated information
P.2.
The grain is case bonded. If not, the entire HIL would be exposed to the high-temperature gas from the beginning of burning, and the simulation would be unnecessary
P.3.
There is no hole or flame-retardant material inside the grain
P.4.
The thickness of the HIL on the same radial section is uniform. Otherwise, manufacturing of the HIL and the grain would become difficult
P.5.
The shape change of the HIL caused by erosion is negligible compared to its overall scaleAccording to the associated information, we define the set of all triangles on the initial burning surface asB, and the set of the triangles on the inner surface of HIL as H.LettingTEx→ be the duration before the flame arriving point x→, we have
(1)Tx→=∫Px→1rt,s→ds,where Px→ is the shortest flame traveling path from B to x→ and rt,s→ is the transient burning rate at the burning location on Px→. Since the burning rate has spatial dependency, the shortest path Px→ may not be a straight line in a general case, where nonuniform burnback models have to be integrated. For the SRMs which can be approximately handled by the uniform burnback model, Px→ is always the shortest line segment connecting B and x→, so the solution of equation (1) can be greatly simplified (see the discussion in Section 2.2).The burning ratert,x→ is decided by the widely known exponential burning rate formula [15]:
(2)rt,x→=Ax→pt,x→nx→,where Ax→ and nx→ are decided by the local propellant material property and initial temperature and pt,x→ is the current gas pressure near x→. Usually, the values of Ax→ and nx→ can be determined by burning rate experiments and are constant among different SRMS, while pt,x→ has to be computed by steady or transient inner ballistic models.For the HIL, the ablation process starts the moment the HIL is exposed to the gas. Lettingte be the elapsed time after exposure, the erosion model at x→ can be summarized by the following time-dependent function:
(3)RRte=MEHtcdtte,Htadtte,Mcte,HM,te,where RRte denotes the recession rate at time te, Htcdt and Htadt denote the heat transfer caused by conduction and advection, respectively, Mc denotes the mechanical effect of condensed particles, and HM stands for the material property of the HIL. The last parameter te implies that function ME has a time-memory property, i.e., ME may output different values at different moments even if Htcdt, Htadt, and Mc remain identical.LetTend be the moment when the SRM stops working. The thickness of the eroded HIL can be computed as follows:
(4)Ex→=∫0Tend−TEx→RRtedte.In order to prevent the overheating of SRM cases, the thickness of the residual HIL atTend must not be too small. The ideal residual thickness can be decided according to the material properties of the SRM. Letting L0 be the minimum allowed residual thickness, we have the following to calculate the optimal thickness of HIL:
(5)Lx→=L0+Ex→.


## 2.2. Model Specialization

Section2.1 presents the universal model which theoretically applies to any SRM. Practically, the performance of most SRMs can be approximately simulated by the uniform burning rate model [9, 10], which assumes that the change in pt,x→ with the position is negligible and thus uses a consistent transient burning rate across the entire grain.By removing the spatial dependency from the burning rate formula, equation (1) can be simplified as the following form:
(6)Dx→,B=∫0Tx→Aptndt,where Dx→,B is the length of the shortest line segment connecting B to x→. There have been some high-performance minimum distance algorithms, including the Havoc algorithm [16] employed in [10], the DiFi algorithm [17] employed in [12], and the Axis-Aligned Bounding Box (AABB) tree method [18]. Both Havoc and DiFi rely on the Voronoi diagram [19] to query the minimum distance. In this article, only the points on H are of our concern, so there are fewer points to process compared to full burnback simulations [10, 12]. Our tests show that the AABB tree method outperforms the Voronoi-type methods when the number of points to be queried is relatively small.


## 3. Numerical Implementation

Here, we present our implementation to discretize and solve the above-established designing model numerically. The solution process can be divided into the following steps:(1)
Discretize the input and find all points for whichTEx→ has to be computed
(2)
Compute values ofTEx→ by the uniform or nonuniform burnback model
(3)
For each of the points, compute the required HIL thickness by the chosen erosion modelIn subsections3.1, 3.2, and 3.3, we provide detailed discussion about the above steps. Since the uniform and the nonuniform burnback models require different operations, their implementations are discussed separately.


### 3.1. Model Discretization

According to the above-discussed property P.4, for each radial section one single value ofL is enough to guide the design. We thus discretize the HIL into ns radial sections, and take nr sample points from each of the sections. ns and nr are large enough to cover any possible geometric features. The discretizing and sampling process can be summarized by Figure 2, where the blue rectangles denote the sections to compute. The lower part of Figure 2 demonstrates a radial section, and the red nodes mark the sampling points.Figure 2
Discretizing and sampling.For the uniform burnback model, which uses an ordinary explicit geometric expression, the coordinates of the sample points can be computed via analytical geometric methods, and the details will be discussed in Section3.2.1. Since the level set method, which is used as the burnback model in this article, does not directly handle explicit geometric elements like sample points, we have developed a “probe vertex” mechanism to imitate the sample points, and the details will be presented in Section 3.3.The above section-by-section discretizing scheme allows us to further simplify the solution process. For each sectionSi, the ideal HIL thickness equals the maximum L value obtained from all its belonging sample points. Letting the set of all sampling points taken from Si be Ssi, the minimum allowed thickness on Si can be calculated as follows:
(7)LSi=max∀x→∈SsiLx→=L0+max∀x→∈SsiEx→.SinceRRte is a nonnegative function, Ex→ monotonically decreases as TEx→ grows. Consequently, equation (7) can be further simplified:
(8)LSi=L0+Eargmin∀x→∈SsiTEx→.Compared to equation (7), equation (8) computes Ex→for only one sample point, evading the numerical integral operations of all other sample points. The general flow of the discretization and the ideal thickness reduction operation can be summarized by Figure 3.Figure 3
Overall flow.


### 3.2. Uniform Burnback Model

#### 3.2.1. Sampling

LettingXminH,XmaxH be the x coordinate range of H, the equation of the ith radial section Si is x=XminH+i−1Δx, where Δx=XmaxH−XminH/ns−1. The intersection of Si and H can be computed by various algorithms, including the AABB tree intersection code provided by [18]. However, there are a few properties that can be utilized to further simplify the computation, which are unlikely to be made use of in general-purpose implementations:
(1)
EachSi is an infinite plane perpendicular to the x axis. The y and z items can therefore be eliminated from the plane equation. In this situation, the intersection of any line and Si can be located via three independent univariate interpolations with respect to x, thus evading the 3 × 3 matrix operation in general approaches
(2)
The facets inH generally distribute perpendicular to Si, so a large portion of the facets can be excluded from the intersecting computation according to their x coordinates. The AABB tree approach [18] contains a similar operation, but the check is unnecessarily performed on all three axes
(3)
Since the sampling is performed on radial sections, all facets parallel toSi can be safely ignored. In general methods, there have to be some extra operations to handle this caseMaking full use of the above properties, we have built a parallelized code, which operates much faster compared to the AABB tree intersection code [18]. The chosen parallel computing framework is the Intel TBB [20], a light-weight framework to distribute tasks on multicore CPUs. The complete process of our code is demonstrated in Procedure 1. In steps 7 and 9 of the Procedure, most nonintersecting triangles to Si are naturally excluded. In step 11, looking up an identifier variable in a table is generally believed to be faster than a series of if-else operations. The two core parts of the Procedure, namely building the intersection polyline and locating the sampling points, are well parallelized. The performance comparison of our code to the AABB tree version is provided in Section 4.


#### 3.2.2. Minimum Distance Querying

As described in Section2, we use the AABB tree distance querying code provided by [18] to query the minimum distance from B to the sampling points. Since the operation is read-only, parallel technology can again be utilized to accelerate the process. The complete flow is demonstrated in Procedure 2.


#### 3.2.3. Exposure Duration Computing

From the above operations, the values ofDx→,B for each sample point are obtained. In order to solve equation (6), the interior ballistic data of the SRM, i.e., pt is still required. The interior ballistic can be computed by any steady or transient zero-dimensional SRM interior ballistic model. In this article, we have used the model developed in [10] to compute the inner ballistic. Once both Dx→,B and pt are ready, Procedure 3 is launched to calculate the exposure duration for each sample point.


### 3.3. Nonuniform Burnback Model

As mentioned above, the level set method is used as the nonuniform burnback model. Details about the level set evolving operations can be found in [12] and are not repeated here. However, the level set evolution process implemented by [12] is performed with a fixed propagation distance step. In order to work with our HIL designing model, the simulation has to be performed on the basis of the time step instead of propagation distance step. The original implementation can be transformed into a time-based one via the following simple steps:
(1)
Using the current gas generation rate data, compute the current inner pressure distribution via any SRM interior ballistic model
(2)
Compute the burning rate distribution on the base of the inner pressure distribution
(3)
Using the intended time step, compute the propagation distance for each mesh vertex
(4)
Perform level set evolution, compute the gas generation rate to be used during the next iterationProcedure 1: Parallel sampling point generation.
Input: Triangle set H, axial discretizing resolution ns, circumferential discretizing resolution nr
1:∀Ti∈H, compute its x coordinate range XminH,XmaxH
2: Sort all triangles inH according to Xmaxi
3: Initialize resulting sampling point setSs=∅
4:Δx=XmaxH−XminH/ns−1
5:for all plane Si at x coordinate Sxi=XminH+i−1Δx serial do
6: Initialize resulting line-segment setRL=∅
7: RemoveTi∣Xmaxi<Sxi from H
8:for all Ti∈H parallel do
9:if Xmini<Sxi then
10:   For each vertexvjj=1,2,3 of Ti, compute its location identifier Li=∑jvxj<Sxi⋅2j
11:   Look upLi in a pre-set table to find the intersecting edges of Ti with Si
12:if There are exactly 2 intersecting edges then
13:    Locate the 2 intersections via univariate interpolation with respect tox
14:    Sort the 2 resulting points by their polar-coordinateθ
15:    Store the line-segment intoRL
16:end if
17:end if
18:end for
19: Sort all line-segments inRL according to the θ polar-coordinate
20: Compute theθ polar-coordinate range θmini,θmaxi of all line-segments in RL
21: OnSi, generate nr rays which start from the center and distribute within θmini,θmaxi
22:for all Lsi∈RL parallel do
23:  Fetch intersecting rays ofLsi according to θ
24:  Locate the sampling points by solving the 2D intersection equations
25:  Store the sampling points intoSS
26:end for
27:end for
28:return SS.Procedure 2: Parallel distance querying.
Input: Triangle set B, sampling point set SS
1: Build a AABB tree object fromB
2:for all sampling sections Si serial do
3: Fetch all sampling points onSi from SS
4: Initialize resulting arrayD
5:for all fetched sampling points parallel do
6:  Query the minimum distance to the AABB tree object, store the result intoD
7:end for
8: Perform parallel reduction to find the minimum distance fromD
9:end for
10:return all computed minimum distances.Procedure 3: Compute exposure duration.
Input:Dx→,B for each sample point, pt data of the SRM
1: Discretizept into a sequence of length np such that ∀i<np,Pi=pti
2: Build the propagation distance data sequenceDp of length np
3:∀i<np, assign Dp∣1=t1rP1, Dp∣i=Dp∣i−1+tirPi
4:for all sampling sections Si parallel do
5: In sequenceDp, find the first i such that Dp∣i−1<DSi,B, Dp∣i≥DSi,B
6: Perform a linear interpolation within intervalti−1,ti, find TESi such that DpTESi=DSi,B,
7:end for
8:return For all the sample points, Tend−TESi is the exposure duration.The above enhancements enable the level set burnback model to calculate the shape of the flame front at the end of each time interval. In order to finally produce the exposure duration data we need, the moment when the flame front bypasses each point on the HIL has to be recorded.Since the level set method expressesB and H by field functions instead of triangles, the above-discussed sampling point mechanism does not apply here, and we introduce “probe vertex” as a replacement. Among the vertices in the level set grid, those inside the HIL are defined as probes. As the flame expands and bypasses one of the probe vertices, the field function on this probe would switch its sign. The moment is recorded as the TEx→ of this probe point. Considering that the main axis of the SRM is aligned with one of the grid axes, the TE values fetched from a slice of the grid are equivalents of the TE values fetched from a sampling section. Like uniform models, the minimum one among the TE values fetched from a slice is taken as the final output of the nonuniform burnback model.


## 3.1. Model Discretization

According to the above-discussed property P.4, for each radial section one single value ofL is enough to guide the design. We thus discretize the HIL into ns radial sections, and take nr sample points from each of the sections. ns and nr are large enough to cover any possible geometric features. The discretizing and sampling process can be summarized by Figure 2, where the blue rectangles denote the sections to compute. The lower part of Figure 2 demonstrates a radial section, and the red nodes mark the sampling points.Figure 2
Discretizing and sampling.For the uniform burnback model, which uses an ordinary explicit geometric expression, the coordinates of the sample points can be computed via analytical geometric methods, and the details will be discussed in Section3.2.1. Since the level set method, which is used as the burnback model in this article, does not directly handle explicit geometric elements like sample points, we have developed a “probe vertex” mechanism to imitate the sample points, and the details will be presented in Section 3.3.The above section-by-section discretizing scheme allows us to further simplify the solution process. For each sectionSi, the ideal HIL thickness equals the maximum L value obtained from all its belonging sample points. Letting the set of all sampling points taken from Si be Ssi, the minimum allowed thickness on Si can be calculated as follows:
(7)LSi=max∀x→∈SsiLx→=L0+max∀x→∈SsiEx→.SinceRRte is a nonnegative function, Ex→ monotonically decreases as TEx→ grows. Consequently, equation (7) can be further simplified:
(8)LSi=L0+Eargmin∀x→∈SsiTEx→.Compared to equation (7), equation (8) computes Ex→for only one sample point, evading the numerical integral operations of all other sample points. The general flow of the discretization and the ideal thickness reduction operation can be summarized by Figure 3.Figure 3
Overall flow.


## 3.2. Uniform Burnback Model

### 3.2.1. Sampling

LettingXminH,XmaxH be the x coordinate range of H, the equation of the ith radial section Si is x=XminH+i−1Δx, where Δx=XmaxH−XminH/ns−1. The intersection of Si and H can be computed by various algorithms, including the AABB tree intersection code provided by [18]. However, there are a few properties that can be utilized to further simplify the computation, which are unlikely to be made use of in general-purpose implementations:
(1)
EachSi is an infinite plane perpendicular to the x axis. The y and z items can therefore be eliminated from the plane equation. In this situation, the intersection of any line and Si can be located via three independent univariate interpolations with respect to x, thus evading the 3 × 3 matrix operation in general approaches
(2)
The facets inH generally distribute perpendicular to Si, so a large portion of the facets can be excluded from the intersecting computation according to their x coordinates. The AABB tree approach [18] contains a similar operation, but the check is unnecessarily performed on all three axes
(3)
Since the sampling is performed on radial sections, all facets parallel toSi can be safely ignored. In general methods, there have to be some extra operations to handle this caseMaking full use of the above properties, we have built a parallelized code, which operates much faster compared to the AABB tree intersection code [18]. The chosen parallel computing framework is the Intel TBB [20], a light-weight framework to distribute tasks on multicore CPUs. The complete process of our code is demonstrated in Procedure 1. In steps 7 and 9 of the Procedure, most nonintersecting triangles to Si are naturally excluded. In step 11, looking up an identifier variable in a table is generally believed to be faster than a series of if-else operations. The two core parts of the Procedure, namely building the intersection polyline and locating the sampling points, are well parallelized. The performance comparison of our code to the AABB tree version is provided in Section 4.


### 3.2.2. Minimum Distance Querying

As described in Section2, we use the AABB tree distance querying code provided by [18] to query the minimum distance from B to the sampling points. Since the operation is read-only, parallel technology can again be utilized to accelerate the process. The complete flow is demonstrated in Procedure 2.


### 3.2.3. Exposure Duration Computing

From the above operations, the values ofDx→,B for each sample point are obtained. In order to solve equation (6), the interior ballistic data of the SRM, i.e., pt is still required. The interior ballistic can be computed by any steady or transient zero-dimensional SRM interior ballistic model. In this article, we have used the model developed in [10] to compute the inner ballistic. Once both Dx→,B and pt are ready, Procedure 3 is launched to calculate the exposure duration for each sample point.


## 3.2.1. Sampling

LettingXminH,XmaxH be the x coordinate range of H, the equation of the ith radial section Si is x=XminH+i−1Δx, where Δx=XmaxH−XminH/ns−1. The intersection of Si and H can be computed by various algorithms, including the AABB tree intersection code provided by [18]. However, there are a few properties that can be utilized to further simplify the computation, which are unlikely to be made use of in general-purpose implementations:
(1)
EachSi is an infinite plane perpendicular to the x axis. The y and z items can therefore be eliminated from the plane equation. In this situation, the intersection of any line and Si can be located via three independent univariate interpolations with respect to x, thus evading the 3 × 3 matrix operation in general approaches
(2)
The facets inH generally distribute perpendicular to Si, so a large portion of the facets can be excluded from the intersecting computation according to their x coordinates. The AABB tree approach [18] contains a similar operation, but the check is unnecessarily performed on all three axes
(3)
Since the sampling is performed on radial sections, all facets parallel toSi can be safely ignored. In general methods, there have to be some extra operations to handle this caseMaking full use of the above properties, we have built a parallelized code, which operates much faster compared to the AABB tree intersection code [18]. The chosen parallel computing framework is the Intel TBB [20], a light-weight framework to distribute tasks on multicore CPUs. The complete process of our code is demonstrated in Procedure 1. In steps 7 and 9 of the Procedure, most nonintersecting triangles to Si are naturally excluded. In step 11, looking up an identifier variable in a table is generally believed to be faster than a series of if-else operations. The two core parts of the Procedure, namely building the intersection polyline and locating the sampling points, are well parallelized. The performance comparison of our code to the AABB tree version is provided in Section 4.


## 3.2.2. Minimum Distance Querying

As described in Section2, we use the AABB tree distance querying code provided by [18] to query the minimum distance from B to the sampling points. Since the operation is read-only, parallel technology can again be utilized to accelerate the process. The complete flow is demonstrated in Procedure 2.


## 3.2.3. Exposure Duration Computing

From the above operations, the values ofDx→,B for each sample point are obtained. In order to solve equation (6), the interior ballistic data of the SRM, i.e., pt is still required. The interior ballistic can be computed by any steady or transient zero-dimensional SRM interior ballistic model. In this article, we have used the model developed in [10] to compute the inner ballistic. Once both Dx→,B and pt are ready, Procedure 3 is launched to calculate the exposure duration for each sample point.


## 3.3. Nonuniform Burnback Model

As mentioned above, the level set method is used as the nonuniform burnback model. Details about the level set evolving operations can be found in [12] and are not repeated here. However, the level set evolution process implemented by [12] is performed with a fixed propagation distance step. In order to work with our HIL designing model, the simulation has to be performed on the basis of the time step instead of propagation distance step. The original implementation can be transformed into a time-based one via the following simple steps:
(1)
Using the current gas generation rate data, compute the current inner pressure distribution via any SRM interior ballistic model
(2)
Compute the burning rate distribution on the base of the inner pressure distribution
(3)
Using the intended time step, compute the propagation distance for each mesh vertex
(4)
Perform level set evolution, compute the gas generation rate to be used during the next iterationProcedure 1: Parallel sampling point generation.
Input: Triangle set H, axial discretizing resolution ns, circumferential discretizing resolution nr
1:∀Ti∈H, compute its x coordinate range XminH,XmaxH
2: Sort all triangles inH according to Xmaxi
3: Initialize resulting sampling point setSs=∅
4:Δx=XmaxH−XminH/ns−1
5:for all plane Si at x coordinate Sxi=XminH+i−1Δx serial do
6: Initialize resulting line-segment setRL=∅
7: RemoveTi∣Xmaxi<Sxi from H
8:for all Ti∈H parallel do
9:if Xmini<Sxi then
10:   For each vertexvjj=1,2,3 of Ti, compute its location identifier Li=∑jvxj<Sxi⋅2j
11:   Look upLi in a pre-set table to find the intersecting edges of Ti with Si
12:if There are exactly 2 intersecting edges then
13:    Locate the 2 intersections via univariate interpolation with respect tox
14:    Sort the 2 resulting points by their polar-coordinateθ
15:    Store the line-segment intoRL
16:end if
17:end if
18:end for
19: Sort all line-segments inRL according to the θ polar-coordinate
20: Compute theθ polar-coordinate range θmini,θmaxi of all line-segments in RL
21: OnSi, generate nr rays which start from the center and distribute within θmini,θmaxi
22:for all Lsi∈RL parallel do
23:  Fetch intersecting rays ofLsi according to θ
24:  Locate the sampling points by solving the 2D intersection equations
25:  Store the sampling points intoSS
26:end for
27:end for
28:return SS.Procedure 2: Parallel distance querying.
Input: Triangle set B, sampling point set SS
1: Build a AABB tree object fromB
2:for all sampling sections Si serial do
3: Fetch all sampling points onSi from SS
4: Initialize resulting arrayD
5:for all fetched sampling points parallel do
6:  Query the minimum distance to the AABB tree object, store the result intoD
7:end for
8: Perform parallel reduction to find the minimum distance fromD
9:end for
10:return all computed minimum distances.Procedure 3: Compute exposure duration.
Input:Dx→,B for each sample point, pt data of the SRM
1: Discretizept into a sequence of length np such that ∀i<np,Pi=pti
2: Build the propagation distance data sequenceDp of length np
3:∀i<np, assign Dp∣1=t1rP1, Dp∣i=Dp∣i−1+tirPi
4:for all sampling sections Si parallel do
5: In sequenceDp, find the first i such that Dp∣i−1<DSi,B, Dp∣i≥DSi,B
6: Perform a linear interpolation within intervalti−1,ti, find TESi such that DpTESi=DSi,B,
7:end for
8:return For all the sample points, Tend−TESi is the exposure duration.The above enhancements enable the level set burnback model to calculate the shape of the flame front at the end of each time interval. In order to finally produce the exposure duration data we need, the moment when the flame front bypasses each point on the HIL has to be recorded.Since the level set method expressesB and H by field functions instead of triangles, the above-discussed sampling point mechanism does not apply here, and we introduce “probe vertex” as a replacement. Among the vertices in the level set grid, those inside the HIL are defined as probes. As the flame expands and bypasses one of the probe vertices, the field function on this probe would switch its sign. The moment is recorded as the TEx→ of this probe point. Considering that the main axis of the SRM is aligned with one of the grid axes, the TE values fetched from a slice of the grid are equivalents of the TE values fetched from a sampling section. Like uniform models, the minimum one among the TE values fetched from a slice is taken as the final output of the nonuniform burnback model.


## 4. Validation and Profiling

Here, we validate the proposed approach and compare the result with that of the traditional method, which uses the constant burning rate presumption.The core idea of this article is to accurately predict the HIL erosion effect and further produce the optimal HIL design in the sense of thickness and weight. The straightforward way to demonstrate the advantage of the proposed designing model is to compare the net weight of the designed HILs with the ones designed by traditional methods. However, the thickness of the HIL is largely affected by the chosen safety factor, and we cannot conclude that one approach is better only because its output produces lighter HIL designing, for the lighter designing may not be safe enough.In this article, we use the following criteria to judge the HIL design method:(1)
The accuracy of HIL erosion prediction. The ideal thickness is computed on the base of the prediction of erosion, so the approach providing higher predicting accuracy is naturally better. The exact HIL erosion data can be obtained by experiments or a full recession simulation performed on both the grain and the HIL
(2)
The unified weight of the designed HIL. In order to eliminate the influence of the safety factor, we scale the thickness of the designs such that the safety factors of both designs are exactly 1. The scaling operation is done as follows:(9)iw=argmin0<i<nsLdSiLeSi,Ld′Si=LdSiLeSiwLdSiw,where LdSi is the thickness function of the designs to be judged, LeSi is the actual minimum required thickness to safely resist heat, and Ld′Si is the unified thickness used to judge the designs. This criterion essentially judges whether some ablation material is “wasted” at positions where the erosion is relatively insignificant.It can be seen that both criteria require exact erosion data. In the discussion below, an extended level set burnback simulation, which treats the HIL as another kind of propellant that its burn rate function equals the erosion rate function, is performed on a sample SRM to provide the exact erosion data. The burning rate distribution is governed as follows:(10)RRt=Ax→ptnx→,x→belongstothegrain,RRcmax0,t−TEx→−Tl,x→belongstotheHIL.After such a level set simulation on the SRM, a field functionϕ containing an implicit expression of the residual HIL is available. Among all the probe vertices within a slice, the minimum ϕ value indicates the thickness of the eroded HIL. The residual thickness can be calculated by subtracting the eroded thickness from the designed thickness. Note that the extended level set simulation is different from the above-proposed nonuniform designing model because performing such an extended simulation requires a much finer grid to capture the geometric characteristics of the HIL.Besides, experiment data on some positions are also provided to validate the extended level set burnback simulation. The sample SRM contains only one kind of propellant and its slenderness ratio is relatively small, therefore the uniform burnback model is integrated to design the HIL. The result of the extended level set simulation and the experiment are provided later in Section4.1.3.For the nonuniform version of the approach, due to the fact that there is no publicly available data and that the evolving accuracy of the employed level set code is already thoroughly verified [12], we only validate the probe vertex mechanism via a simple nonuniform grain having a clear analytical solution. In addition, the computing efficiency of the code is profiled to demonstrate the performance of the developed code.


### 4.1. Uniform Approach

#### 4.1.1. Sample SRM

The sample SRM is demonstrated in Figure4, where the nozzle is mounted to the right side of the figure. The head and the tail part of the grain each contains a group of 10 fins. The radius of the SRM is 300 mm, and the length is 1900 mm (the ellipsoid parts are included). The radius of the inner hole is 85 mm. The height of the fins is 220 mm.Figure 4
The grain of the sample SRM.
(a)
The sectional view
(b)
The initial burn surfaceThe burning rate of the propellant is decided by equation (2). Burning rate experiments at room temperature show that the burning rate constants of the propellant used in the sample SRM are n=0.4, A=5.2015mm/MPan. The radius of the nozzle throat is 50 mm. In Figure 5, we have demonstrated the interior flow field of the sample SRM. It can be seen that the difference of the static pressure across the entire flow field is insignificant (the pressure range is 6.63–6.94 MPa and the relative burning rate variation is within 1.83%), so applying the uniform burnback model on the sample SRM is practical. In Figure 6, the red line demonstrated the inner ballistic acquired via the uniform burning rate assumption and the zero-dimensional inner ballistic model. The blue line plots the fluctuation of the propellant burning rate computed from equation (2).Figure 5
Interior flow field of the sample SRM.Figure 6
Inner ballistic and transient burning rate of the sample SRM.The sample SRM uses nitrile butadiene rubber (NBR) as the insulation material, and the erosion characteristics are obtained from ablation experiments. After being exposed in the stream composed of hot gas and condensed particles, the HIL keeps its original shape within a short duration, and then its boundary starts to recess at an approximately constant speed. This behavior is explained by Yang et al. in [6]: the shape of the HIL remains constant before the carbonaceous char layer appears and starts to move inward. Generally, the erosion model of the insulation material used in the sample SRM can be summarized as follows:
(11)RRte=0,te<Tl,RRc,te≥Tl,where Tl is the duration before the carbonaceous char layer starts to move inward. The constants measured from the experiment data are Tl=1.3s and RRc=0.35mm/s.Substituting equation (11) into equation (4) and equation (5), we have the following compact equation to design the thickness of the HIL:
(12)Lx→=L0+0,Tend−TEx→<Tl,RRcTend−TEx→,Tend−TEx→≥Tl.Other erosion models can be simply integrated by replacing theRRcTend−TEx→ item with corresponding numerical integrations.


#### 4.1.2. Burnback Model

The SRM is then fed to the above-implemented code. In Figure7, the generated sampling points, colored according to their D values, are demonstrated. The D value and exposure duration on each sampling section is plotted in Figure 8 by the blue line and the red line, respectively.Figure 7
D on a sampling point set.Figure 8
D and exposure duration on sampling sections (ns=1024 and nr=180).Besides, a manual CAD measuring operation is accomplished on 29 chosen sections of the grain (the chosen sections are denser within the head and tail parts compared to the center segment). In Figure8, we have also plotted the manually measured D values from the CAD platform. It can be seen that every D data point measured from CAD exactly falls on the curve outputted by the proposed method, but the CAD measured data misses some feature points due to the insufficiency of sampling sections. The constant burning rate used for the CAD-based traditional method is computed by dividing the maximum D value among all sections (215 mm) by the working duration of the SRM (18.3729 s), and the result is 11.7020 mm/s. Using the data, the exposure duration prediction of the CAD-based method is computed and plotted in Figure 8. It can be seen that the traditional approach tends to give smaller predictions of the exposure duration near the ellipsoid segments of the sample SRM. This is due to the fact that the sample grain burns relatively faster in its first half working time, and thus, it exposes some of the HIL earlier than the prediction made on the basis of a constant burning rate. When the pressure exponent n in equation (2) becomes larger, the difference would grow more significant. Due to similar reasons, for the SRMs which burn slower during its front stage, the traditional approach would output a longer prediction.


#### 4.1.3. Thickness Prediction

Next, two HILs are designed via the CAD approach and the proposed approach, respectively. The suggested thickness of both approaches can be directly computed from equation (12) (the required L0 is set as 3 mm according to engineering experience) and the above-acquired exposure duration data. In Figure 9, we have plotted the resulting thickness of both approaches and the actually required thickness (eroded thickness plus L0) data obtained from the abovementioned extended level set simulation. It can be seen that the data acquired from the extended simulation fits well with the prediction of the proposed approach. The only two differences of the “proposed approach” and the “level set prediction” curves lie in the region where the eroded thickness changes dramatically, where the “level set prediction” curve is smoother due to the level set method’s nature of smearing out subcell geometric features. The “CAD approach” curve, on the other hand, suggests thinner thickness in several regions, implying that the approach failed to accurately predict the erosion. The proposed method satisfies the abovementioned criterion (1).Figure 9
ResultingL values on sampling sections.From the tail dome of an actual SRM sample which had experienced a thrust experiment, we have measured the actual eroded thickness of the HIL from 3 points (see Figure10, a thorough measurement is not possible because damaging the dome is not allowed). The acquired data, listed in Table 1, shows that the output of the extended level set simulation is reliable and can be used to replace the experiment data.Figure 10
Positions of the 3 measured points.Table 1
Comparison of measured and predicted erosion.
Point description
Eroded thickness (mm)
Measured
Predicted
A: adjacent to the nozzle, exposed at the start
5.95
5.9720
B: on the front edge of the dome, exposed at 6.24 s
3.78
3.7915
C: on the front edge of the dome, exposed at 8.28 s
2.99
3.0775Using equation (9), we can scale both designing schemes such that the safety factors are exactly 1 with respect to the prediction of the level set simulation. The scaled designs are outlined in Figure 11. It can be seen that the design produced by the CAD approach is generally thicker after the unifying operation, meaning that some of insulating material does not contribute to the overall safety factor. The volumes of the unified designs (calculated via analytical geometric methods) are 1.30e−4m3 and 1.35e−4m3, respectively. The design produced by the proposed approach is 3.85% lighter compared to that of the traditional approach.Figure 11
Unified designs.


#### 4.1.4. Effect of Discretizing Resolutions

The sampling resolution used in the above process isns=1024, nr=180, which produces an accurate D compared to the CAD measured value. In order to characterize the effect of sampling resolution on the error, we have launched more tests using various ns and nr combinations and various SRM models. The results can be summarized by the following conclusions:
(1)
IfH is a revolving surface, nr≥5Nas usually provides enough accuracy, where Nas is the number of axial sweeping geometric features (e.g., fins); setting a larger nr while leaving ns unmoved produces little improvement in accuracy
(2)
IfH is not a revolving surface, nr should be set such that the spacing between neighboring sampling points is smaller than the minimum characteristic scale of H
(3)
ns should be set such that the spacing between neighboring sampling sections is smaller than the minimum characteristic scale among the SRM


### 4.2. Nonuniform Approach

The nonuniform burnback model is validated by comparing the computation output to the analytical result. As discussed above, we only need to check if the probe mechanism works as expected. The SRM used in the verification is demonstrated in Figure12, where the red and blue thick lines mark B and H, respectively. The dotted area is filled by a slower-burning propellant, and the crossed area is filled by a faster-burning propellant. The interface between the two kinds of propellants is exactly located at the center of the SRM. The head and tail part of the SRM grain is ellipsoid.Figure 12
A sample grain to verify the nonuniform burnback model.According to [21] and the ellipse equation, the following analytical result can be derived (assuming l≫a−h):
(13)Xm=b1−h2a2,Rsf=rsrf,Xs=tansin−1Rsf1−Rsfa−h,TEx=0,x<−l+Xm,orx>l+Xm,1rsa1−x+12b2−h,−l+Xm≤x<−l,a−hrs,−l≤x<−Xs,a−h1Xs1rf−1rsx+1rf,−Xs≤x<0,a−hrf,0≤x<l,1rfa1−x−12b2−h,l≤x<l+Xm,where rs and rf are the constant burning rates of the slower- and faster-burning propellants, respectively.Lettingl=100mm, a=40mm, b=15mm, h=15mm, rs=6mm/s, and rf=9mm/s, the outputs of the nonuniform burnback model are collected in Figure 13. Equation (12) is again used as the erosion model. Obviously, the simulation outputs fit well the analytical result, and the error tends to be negligible as the level set mesh is refined.Figure 13
TE obtained from the nonuniform burnback model.


### 4.3. Performance Profiling

There are mainly four acceleration technologies utilized in the proposed designing process:(1)
The efficient sampling point-generating process, which makes use of the special geometric properties of SRMs to speed up the sampling (Section3.2.1)
(2)
The parallel computing technology utilized in Procedure1
(3)
The parallel computing technology utilized in Procedure2
(4)
The parallel computing technology utilized in the nonuniform burnback modelProfiling data of the first three technologies are demonstrated in the discussion below. For the nonuniform burnback model, detailed profiling data is provided in [12].Under normal circumstances, it takes less than 5 s for our optimized code to perform a 4-thread simulation using the uniform burnback model, while the unoptimized serial version costs more than 30 s.


#### 4.3.1. Efficient Sampling Point-Generating Process: Procedure1

Discretizing the CAD model shown by Figure4 using different controlling parameters, three triangular models containing 1336, 3497, and 4635 triangles, respectively, are generated. Then, we generate the sampling points via the AABB tree approach [18] and Procedure 1, respectively. Repeating the computation and averaging the execution time, we get the profiling result shown by Figure 14, where nr=180, ns=128,256,512,1024, and the “FAST” labels mark the data measured from Procedure 1. Here, the parallel computing in Procedure 1 is disabled by setting the maximum allowed number of threads to 1. The profiling is performed on an Intel Core i5 3230M CPU and DDR3 1.33 GHz RAM.Figure 14
Serial sampling point-generating profiling data.Obviously, the efficiency of the proposed algorithm is significantly higher than that of the AABB tree approach, with whatever discretizing resolution. More precisely, Procedure1 saves roughly 35% of the execution time compared to the AABB tree approach.


#### 4.3.2. Parallel Generating of Sampling Points: Procedure1

Here, we measure the improvement introduced by the parallel computing in Procedure1. The triangular model having 4635 triangles is employed as the testing model. Accuracy parameters are set as nr=180 and ns=128,256,512,1024. The performance of the code is measured by the metric “Sampling point-generating rate” plotted in Figure 15, where the unit “Mspps” means “million sampling points per second.” The hardware platform is an Intel Xeon E5 2690 v4 CPU with a DDR4 2.4 GHz RAM.Figure 15
Parallel sampling point-generating profiling data.In Figure15, the numbered legends denote the ns used to generate the curve. The “Ideal” curve shows the theoretical upper performance bound when parallel cost is negligible. The curves show that the generating rates grow as the number of threads rises, but the growing rates are slower than the ideal linear growth. This is an inevitable issue for any parallel code, since some computation power has to be spent on organizing the threads, and some threads may be idle at the end of a parallel loop since there is no more data to be processed. Moreover, the generating rate curve slightly shifts upward as ns grows. This phenomenon is caused by the nearly constant initiating cost in the code: when there are more sampling points to be generated, the negative impact of the initiating cost would naturally be smaller after averaging.


#### 4.3.3. Parallel Distance Querying: Procedure2

Using the same geometric model and hardware platform as Section4.3.2, we have obtained the parallel performance of Procedure 2 shown by Figure 16. Since most of the steps in Procedure 2 are parallelized, there is less initiating cost in Procedure 2 compared to Procedure 1, so the difference caused by different ns values is negligible. Figure 16 shows the averaged distance querying rate under various discretizing resolutions. It can be seen that the querying rate grows nearly linearly as the number of threads increases.Figure 16
Parallel distance querying profiling data.


## 4.1. Uniform Approach

### 4.1.1. Sample SRM

The sample SRM is demonstrated in Figure4, where the nozzle is mounted to the right side of the figure. The head and the tail part of the grain each contains a group of 10 fins. The radius of the SRM is 300 mm, and the length is 1900 mm (the ellipsoid parts are included). The radius of the inner hole is 85 mm. The height of the fins is 220 mm.Figure 4
The grain of the sample SRM.
(a)
The sectional view
(b)
The initial burn surfaceThe burning rate of the propellant is decided by equation (2). Burning rate experiments at room temperature show that the burning rate constants of the propellant used in the sample SRM are n=0.4, A=5.2015mm/MPan. The radius of the nozzle throat is 50 mm. In Figure 5, we have demonstrated the interior flow field of the sample SRM. It can be seen that the difference of the static pressure across the entire flow field is insignificant (the pressure range is 6.63–6.94 MPa and the relative burning rate variation is within 1.83%), so applying the uniform burnback model on the sample SRM is practical. In Figure 6, the red line demonstrated the inner ballistic acquired via the uniform burning rate assumption and the zero-dimensional inner ballistic model. The blue line plots the fluctuation of the propellant burning rate computed from equation (2).Figure 5
Interior flow field of the sample SRM.Figure 6
Inner ballistic and transient burning rate of the sample SRM.The sample SRM uses nitrile butadiene rubber (NBR) as the insulation material, and the erosion characteristics are obtained from ablation experiments. After being exposed in the stream composed of hot gas and condensed particles, the HIL keeps its original shape within a short duration, and then its boundary starts to recess at an approximately constant speed. This behavior is explained by Yang et al. in [6]: the shape of the HIL remains constant before the carbonaceous char layer appears and starts to move inward. Generally, the erosion model of the insulation material used in the sample SRM can be summarized as follows:
(11)RRte=0,te<Tl,RRc,te≥Tl,where Tl is the duration before the carbonaceous char layer starts to move inward. The constants measured from the experiment data are Tl=1.3s and RRc=0.35mm/s.Substituting equation (11) into equation (4) and equation (5), we have the following compact equation to design the thickness of the HIL:
(12)Lx→=L0+0,Tend−TEx→<Tl,RRcTend−TEx→,Tend−TEx→≥Tl.Other erosion models can be simply integrated by replacing theRRcTend−TEx→ item with corresponding numerical integrations.


### 4.1.2. Burnback Model

The SRM is then fed to the above-implemented code. In Figure7, the generated sampling points, colored according to their D values, are demonstrated. The D value and exposure duration on each sampling section is plotted in Figure 8 by the blue line and the red line, respectively.Figure 7
D on a sampling point set.Figure 8
D and exposure duration on sampling sections (ns=1024 and nr=180).Besides, a manual CAD measuring operation is accomplished on 29 chosen sections of the grain (the chosen sections are denser within the head and tail parts compared to the center segment). In Figure8, we have also plotted the manually measured D values from the CAD platform. It can be seen that every D data point measured from CAD exactly falls on the curve outputted by the proposed method, but the CAD measured data misses some feature points due to the insufficiency of sampling sections. The constant burning rate used for the CAD-based traditional method is computed by dividing the maximum D value among all sections (215 mm) by the working duration of the SRM (18.3729 s), and the result is 11.7020 mm/s. Using the data, the exposure duration prediction of the CAD-based method is computed and plotted in Figure 8. It can be seen that the traditional approach tends to give smaller predictions of the exposure duration near the ellipsoid segments of the sample SRM. This is due to the fact that the sample grain burns relatively faster in its first half working time, and thus, it exposes some of the HIL earlier than the prediction made on the basis of a constant burning rate. When the pressure exponent n in equation (2) becomes larger, the difference would grow more significant. Due to similar reasons, for the SRMs which burn slower during its front stage, the traditional approach would output a longer prediction.


### 4.1.3. Thickness Prediction

Next, two HILs are designed via the CAD approach and the proposed approach, respectively. The suggested thickness of both approaches can be directly computed from equation (12) (the required L0 is set as 3 mm according to engineering experience) and the above-acquired exposure duration data. In Figure 9, we have plotted the resulting thickness of both approaches and the actually required thickness (eroded thickness plus L0) data obtained from the abovementioned extended level set simulation. It can be seen that the data acquired from the extended simulation fits well with the prediction of the proposed approach. The only two differences of the “proposed approach” and the “level set prediction” curves lie in the region where the eroded thickness changes dramatically, where the “level set prediction” curve is smoother due to the level set method’s nature of smearing out subcell geometric features. The “CAD approach” curve, on the other hand, suggests thinner thickness in several regions, implying that the approach failed to accurately predict the erosion. The proposed method satisfies the abovementioned criterion (1).Figure 9
ResultingL values on sampling sections.From the tail dome of an actual SRM sample which had experienced a thrust experiment, we have measured the actual eroded thickness of the HIL from 3 points (see Figure10, a thorough measurement is not possible because damaging the dome is not allowed). The acquired data, listed in Table 1, shows that the output of the extended level set simulation is reliable and can be used to replace the experiment data.Figure 10
Positions of the 3 measured points.Table 1
Comparison of measured and predicted erosion.
Point description
Eroded thickness (mm)
Measured
Predicted
A: adjacent to the nozzle, exposed at the start
5.95
5.9720
B: on the front edge of the dome, exposed at 6.24 s
3.78
3.7915
C: on the front edge of the dome, exposed at 8.28 s
2.99
3.0775Using equation (9), we can scale both designing schemes such that the safety factors are exactly 1 with respect to the prediction of the level set simulation. The scaled designs are outlined in Figure 11. It can be seen that the design produced by the CAD approach is generally thicker after the unifying operation, meaning that some of insulating material does not contribute to the overall safety factor. The volumes of the unified designs (calculated via analytical geometric methods) are 1.30e−4m3 and 1.35e−4m3, respectively. The design produced by the proposed approach is 3.85% lighter compared to that of the traditional approach.Figure 11
Unified designs.


### 4.1.4. Effect of Discretizing Resolutions

The sampling resolution used in the above process isns=1024, nr=180, which produces an accurate D compared to the CAD measured value. In order to characterize the effect of sampling resolution on the error, we have launched more tests using various ns and nr combinations and various SRM models. The results can be summarized by the following conclusions:
(1)
IfH is a revolving surface, nr≥5Nas usually provides enough accuracy, where Nas is the number of axial sweeping geometric features (e.g., fins); setting a larger nr while leaving ns unmoved produces little improvement in accuracy
(2)
IfH is not a revolving surface, nr should be set such that the spacing between neighboring sampling points is smaller than the minimum characteristic scale of H
(3)
ns should be set such that the spacing between neighboring sampling sections is smaller than the minimum characteristic scale among the SRM


## 4.1.1. Sample SRM

The sample SRM is demonstrated in Figure4, where the nozzle is mounted to the right side of the figure. The head and the tail part of the grain each contains a group of 10 fins. The radius of the SRM is 300 mm, and the length is 1900 mm (the ellipsoid parts are included). The radius of the inner hole is 85 mm. The height of the fins is 220 mm.Figure 4
The grain of the sample SRM.
(a)
The sectional view
(b)
The initial burn surfaceThe burning rate of the propellant is decided by equation (2). Burning rate experiments at room temperature show that the burning rate constants of the propellant used in the sample SRM are n=0.4, A=5.2015mm/MPan. The radius of the nozzle throat is 50 mm. In Figure 5, we have demonstrated the interior flow field of the sample SRM. It can be seen that the difference of the static pressure across the entire flow field is insignificant (the pressure range is 6.63–6.94 MPa and the relative burning rate variation is within 1.83%), so applying the uniform burnback model on the sample SRM is practical. In Figure 6, the red line demonstrated the inner ballistic acquired via the uniform burning rate assumption and the zero-dimensional inner ballistic model. The blue line plots the fluctuation of the propellant burning rate computed from equation (2).Figure 5
Interior flow field of the sample SRM.Figure 6
Inner ballistic and transient burning rate of the sample SRM.The sample SRM uses nitrile butadiene rubber (NBR) as the insulation material, and the erosion characteristics are obtained from ablation experiments. After being exposed in the stream composed of hot gas and condensed particles, the HIL keeps its original shape within a short duration, and then its boundary starts to recess at an approximately constant speed. This behavior is explained by Yang et al. in [6]: the shape of the HIL remains constant before the carbonaceous char layer appears and starts to move inward. Generally, the erosion model of the insulation material used in the sample SRM can be summarized as follows:
(11)RRte=0,te<Tl,RRc,te≥Tl,where Tl is the duration before the carbonaceous char layer starts to move inward. The constants measured from the experiment data are Tl=1.3s and RRc=0.35mm/s.Substituting equation (11) into equation (4) and equation (5), we have the following compact equation to design the thickness of the HIL:
(12)Lx→=L0+0,Tend−TEx→<Tl,RRcTend−TEx→,Tend−TEx→≥Tl.Other erosion models can be simply integrated by replacing theRRcTend−TEx→ item with corresponding numerical integrations.


## 4.1.2. Burnback Model

The SRM is then fed to the above-implemented code. In Figure7, the generated sampling points, colored according to their D values, are demonstrated. The D value and exposure duration on each sampling section is plotted in Figure 8 by the blue line and the red line, respectively.Figure 7
D on a sampling point set.Figure 8
D and exposure duration on sampling sections (ns=1024 and nr=180).Besides, a manual CAD measuring operation is accomplished on 29 chosen sections of the grain (the chosen sections are denser within the head and tail parts compared to the center segment). In Figure8, we have also plotted the manually measured D values from the CAD platform. It can be seen that every D data point measured from CAD exactly falls on the curve outputted by the proposed method, but the CAD measured data misses some feature points due to the insufficiency of sampling sections. The constant burning rate used for the CAD-based traditional method is computed by dividing the maximum D value among all sections (215 mm) by the working duration of the SRM (18.3729 s), and the result is 11.7020 mm/s. Using the data, the exposure duration prediction of the CAD-based method is computed and plotted in Figure 8. It can be seen that the traditional approach tends to give smaller predictions of the exposure duration near the ellipsoid segments of the sample SRM. This is due to the fact that the sample grain burns relatively faster in its first half working time, and thus, it exposes some of the HIL earlier than the prediction made on the basis of a constant burning rate. When the pressure exponent n in equation (2) becomes larger, the difference would grow more significant. Due to similar reasons, for the SRMs which burn slower during its front stage, the traditional approach would output a longer prediction.


## 4.1.3. Thickness Prediction

Next, two HILs are designed via the CAD approach and the proposed approach, respectively. The suggested thickness of both approaches can be directly computed from equation (12) (the required L0 is set as 3 mm according to engineering experience) and the above-acquired exposure duration data. In Figure 9, we have plotted the resulting thickness of both approaches and the actually required thickness (eroded thickness plus L0) data obtained from the abovementioned extended level set simulation. It can be seen that the data acquired from the extended simulation fits well with the prediction of the proposed approach. The only two differences of the “proposed approach” and the “level set prediction” curves lie in the region where the eroded thickness changes dramatically, where the “level set prediction” curve is smoother due to the level set method’s nature of smearing out subcell geometric features. The “CAD approach” curve, on the other hand, suggests thinner thickness in several regions, implying that the approach failed to accurately predict the erosion. The proposed method satisfies the abovementioned criterion (1).Figure 9
ResultingL values on sampling sections.From the tail dome of an actual SRM sample which had experienced a thrust experiment, we have measured the actual eroded thickness of the HIL from 3 points (see Figure10, a thorough measurement is not possible because damaging the dome is not allowed). The acquired data, listed in Table 1, shows that the output of the extended level set simulation is reliable and can be used to replace the experiment data.Figure 10
Positions of the 3 measured points.Table 1
Comparison of measured and predicted erosion.
Point description
Eroded thickness (mm)
Measured
Predicted
A: adjacent to the nozzle, exposed at the start
5.95
5.9720
B: on the front edge of the dome, exposed at 6.24 s
3.78
3.7915
C: on the front edge of the dome, exposed at 8.28 s
2.99
3.0775Using equation (9), we can scale both designing schemes such that the safety factors are exactly 1 with respect to the prediction of the level set simulation. The scaled designs are outlined in Figure 11. It can be seen that the design produced by the CAD approach is generally thicker after the unifying operation, meaning that some of insulating material does not contribute to the overall safety factor. The volumes of the unified designs (calculated via analytical geometric methods) are 1.30e−4m3 and 1.35e−4m3, respectively. The design produced by the proposed approach is 3.85% lighter compared to that of the traditional approach.Figure 11
Unified designs.


## 4.1.4. Effect of Discretizing Resolutions

The sampling resolution used in the above process isns=1024, nr=180, which produces an accurate D compared to the CAD measured value. In order to characterize the effect of sampling resolution on the error, we have launched more tests using various ns and nr combinations and various SRM models. The results can be summarized by the following conclusions:
(1)
IfH is a revolving surface, nr≥5Nas usually provides enough accuracy, where Nas is the number of axial sweeping geometric features (e.g., fins); setting a larger nr while leaving ns unmoved produces little improvement in accuracy
(2)
IfH is not a revolving surface, nr should be set such that the spacing between neighboring sampling points is smaller than the minimum characteristic scale of H
(3)
ns should be set such that the spacing between neighboring sampling sections is smaller than the minimum characteristic scale among the SRM


## 4.2. Nonuniform Approach

The nonuniform burnback model is validated by comparing the computation output to the analytical result. As discussed above, we only need to check if the probe mechanism works as expected. The SRM used in the verification is demonstrated in Figure12, where the red and blue thick lines mark B and H, respectively. The dotted area is filled by a slower-burning propellant, and the crossed area is filled by a faster-burning propellant. The interface between the two kinds of propellants is exactly located at the center of the SRM. The head and tail part of the SRM grain is ellipsoid.Figure 12
A sample grain to verify the nonuniform burnback model.According to [21] and the ellipse equation, the following analytical result can be derived (assuming l≫a−h):
(13)Xm=b1−h2a2,Rsf=rsrf,Xs=tansin−1Rsf1−Rsfa−h,TEx=0,x<−l+Xm,orx>l+Xm,1rsa1−x+12b2−h,−l+Xm≤x<−l,a−hrs,−l≤x<−Xs,a−h1Xs1rf−1rsx+1rf,−Xs≤x<0,a−hrf,0≤x<l,1rfa1−x−12b2−h,l≤x<l+Xm,where rs and rf are the constant burning rates of the slower- and faster-burning propellants, respectively.Lettingl=100mm, a=40mm, b=15mm, h=15mm, rs=6mm/s, and rf=9mm/s, the outputs of the nonuniform burnback model are collected in Figure 13. Equation (12) is again used as the erosion model. Obviously, the simulation outputs fit well the analytical result, and the error tends to be negligible as the level set mesh is refined.Figure 13
TE obtained from the nonuniform burnback model.


## 4.3. Performance Profiling

There are mainly four acceleration technologies utilized in the proposed designing process:(1)
The efficient sampling point-generating process, which makes use of the special geometric properties of SRMs to speed up the sampling (Section3.2.1)
(2)
The parallel computing technology utilized in Procedure1
(3)
The parallel computing technology utilized in Procedure2
(4)
The parallel computing technology utilized in the nonuniform burnback modelProfiling data of the first three technologies are demonstrated in the discussion below. For the nonuniform burnback model, detailed profiling data is provided in [12].Under normal circumstances, it takes less than 5 s for our optimized code to perform a 4-thread simulation using the uniform burnback model, while the unoptimized serial version costs more than 30 s.


### 4.3.1. Efficient Sampling Point-Generating Process: Procedure1

Discretizing the CAD model shown by Figure4 using different controlling parameters, three triangular models containing 1336, 3497, and 4635 triangles, respectively, are generated. Then, we generate the sampling points via the AABB tree approach [18] and Procedure 1, respectively. Repeating the computation and averaging the execution time, we get the profiling result shown by Figure 14, where nr=180, ns=128,256,512,1024, and the “FAST” labels mark the data measured from Procedure 1. Here, the parallel computing in Procedure 1 is disabled by setting the maximum allowed number of threads to 1. The profiling is performed on an Intel Core i5 3230M CPU and DDR3 1.33 GHz RAM.Figure 14
Serial sampling point-generating profiling data.Obviously, the efficiency of the proposed algorithm is significantly higher than that of the AABB tree approach, with whatever discretizing resolution. More precisely, Procedure1 saves roughly 35% of the execution time compared to the AABB tree approach.


### 4.3.2. Parallel Generating of Sampling Points: Procedure1

Here, we measure the improvement introduced by the parallel computing in Procedure1. The triangular model having 4635 triangles is employed as the testing model. Accuracy parameters are set as nr=180 and ns=128,256,512,1024. The performance of the code is measured by the metric “Sampling point-generating rate” plotted in Figure 15, where the unit “Mspps” means “million sampling points per second.” The hardware platform is an Intel Xeon E5 2690 v4 CPU with a DDR4 2.4 GHz RAM.Figure 15
Parallel sampling point-generating profiling data.In Figure15, the numbered legends denote the ns used to generate the curve. The “Ideal” curve shows the theoretical upper performance bound when parallel cost is negligible. The curves show that the generating rates grow as the number of threads rises, but the growing rates are slower than the ideal linear growth. This is an inevitable issue for any parallel code, since some computation power has to be spent on organizing the threads, and some threads may be idle at the end of a parallel loop since there is no more data to be processed. Moreover, the generating rate curve slightly shifts upward as ns grows. This phenomenon is caused by the nearly constant initiating cost in the code: when there are more sampling points to be generated, the negative impact of the initiating cost would naturally be smaller after averaging.


### 4.3.3. Parallel Distance Querying: Procedure2

Using the same geometric model and hardware platform as Section4.3.2, we have obtained the parallel performance of Procedure 2 shown by Figure 16. Since most of the steps in Procedure 2 are parallelized, there is less initiating cost in Procedure 2 compared to Procedure 1, so the difference caused by different ns values is negligible. Figure 16 shows the averaged distance querying rate under various discretizing resolutions. It can be seen that the querying rate grows nearly linearly as the number of threads increases.Figure 16
Parallel distance querying profiling data.


## 4.3.1. Efficient Sampling Point-Generating Process: Procedure1

Discretizing the CAD model shown by Figure4 using different controlling parameters, three triangular models containing 1336, 3497, and 4635 triangles, respectively, are generated. Then, we generate the sampling points via the AABB tree approach [18] and Procedure 1, respectively. Repeating the computation and averaging the execution time, we get the profiling result shown by Figure 14, where nr=180, ns=128,256,512,1024, and the “FAST” labels mark the data measured from Procedure 1. Here, the parallel computing in Procedure 1 is disabled by setting the maximum allowed number of threads to 1. The profiling is performed on an Intel Core i5 3230M CPU and DDR3 1.33 GHz RAM.Figure 14
Serial sampling point-generating profiling data.Obviously, the efficiency of the proposed algorithm is significantly higher than that of the AABB tree approach, with whatever discretizing resolution. More precisely, Procedure1 saves roughly 35% of the execution time compared to the AABB tree approach.


## 4.3.2. Parallel Generating of Sampling Points: Procedure1

Here, we measure the improvement introduced by the parallel computing in Procedure1. The triangular model having 4635 triangles is employed as the testing model. Accuracy parameters are set as nr=180 and ns=128,256,512,1024. The performance of the code is measured by the metric “Sampling point-generating rate” plotted in Figure 15, where the unit “Mspps” means “million sampling points per second.” The hardware platform is an Intel Xeon E5 2690 v4 CPU with a DDR4 2.4 GHz RAM.Figure 15
Parallel sampling point-generating profiling data.In Figure15, the numbered legends denote the ns used to generate the curve. The “Ideal” curve shows the theoretical upper performance bound when parallel cost is negligible. The curves show that the generating rates grow as the number of threads rises, but the growing rates are slower than the ideal linear growth. This is an inevitable issue for any parallel code, since some computation power has to be spent on organizing the threads, and some threads may be idle at the end of a parallel loop since there is no more data to be processed. Moreover, the generating rate curve slightly shifts upward as ns grows. This phenomenon is caused by the nearly constant initiating cost in the code: when there are more sampling points to be generated, the negative impact of the initiating cost would naturally be smaller after averaging.


## 4.3.3. Parallel Distance Querying: Procedure2

Using the same geometric model and hardware platform as Section4.3.2, we have obtained the parallel performance of Procedure 2 shown by Figure 16. Since most of the steps in Procedure 2 are parallelized, there is less initiating cost in Procedure 2 compared to Procedure 1, so the difference caused by different ns values is negligible. Figure 16 shows the averaged distance querying rate under various discretizing resolutions. It can be seen that the querying rate grows nearly linearly as the number of threads increases.Figure 16
Parallel distance querying profiling data.


## 5. Conclusions

In this article, we have established a set of numerical models to design HILs. The proposed models are compatible with SRMs having any shape, any dynamic burning rate distribution, and any manner of erosion. Validation results show that the model produces optimal designs for the HILs, saving roughly 3.85% of the mass for HILs. The proposed approach does not rely on CAD platforms and eliminates all manual operations.We have also provided details about the high-performance code implementation of the numerical model. The provided implementation is well parallelized to make full use of modern multicore CPUs. An innovative sampling point-generating algorithm is included in the code, saving 35% of the execution time compared to its open source equivalent. Normal simulations cost less than 5 s on modern 4-thread CPUs.Due to manufacturing process considerations, the proposed designing routine assumes the thickness of HILs to be uniform on the same radial section. However, the sampling point mechanism in our approach can be easily extended to allow outputting the suggested thickness for each point of the HIL. Once such three-dimensional designs become practical, the proposed method can be extended to design even thinner and lighter HILs while ensuring safety.



---

*Source: 5789430-2019-01-21.xml*