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

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.


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 Figure 1. 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.
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][7][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][12][13][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 Section 2 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.

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 scale According to the associated information, we define the set of all triangles on the initial burning surface as B, and the set of the triangles on the inner surface of HIL as H.
Letting T E x be the duration before the flame arriving point x , we have where P x is the shortest flame traveling path from B to x and r t, s is the transient burning rate at the burning   (1) can be greatly simplified (see the discussion in Section 2.2). The burning rate r t, x is decided by the widely known exponential burning rate formula [15]: where A x and n x are decided by the local propellant material property and initial temperature and p t, x is the current gas pressure near x . Usually, the values of A x and n x can be determined by burning rate experiments and are constant among different SRMS, while p t, 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. Letting t e be the elapsed time after exposure, the erosion model at x can be summarized by the following time-dependent function: In order to prevent the overheating of SRM cases, the thickness of the residual HIL at T end must not be too small. The ideal residual thickness can be decided according to the material properties of the SRM. Letting L 0 be the minimum allowed residual thickness, we have the following to calculate the optimal thickness of HIL: 2.2. Model Specialization. Section 2.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 p t, 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: where D x , B is the length of the shortest line segment connecting B to x . There have been some highperformance 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.

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 which T E x has to be computed In subsections 3.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 of L is enough to guide the design. We thus discretize the HIL into n s radial sections, and take n r sample points from each of the sections. n s and n r 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.
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 Section 3.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.

International Journal of Aerospace Engineering
The above section-by-section discretizing scheme allows us to further simplify the solution process. For each section S i , 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 S i be S i s , the minimum allowed thickness on S i can be calculated as follows: Since R R t e is a nonnegative function, E x monotonically decreases as T E x grows. Consequently, equation (7) can be further simplified: Compared to equation (7), equation (8) computes E x 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.

Uniform Burnback Model
The intersection of S i 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) Each S i 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 S i can be located via three independent univariate interpolations with respect to x, thus evading the 3 × 3 matrix operation in general approaches Making 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 S i 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 Section 2, 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.

Exposure Duration
Computing. From the above operations, the values of D x , B for each sample point are obtained. In order to solve equation (6), the interior ballistic data of the SRM, i.e., p t 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 D x , B and p t 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)  if X i min < S i x then 10: For each vertex v j j = 1, 2, 3 of T i , compute its location identifier Look up L i in a pre-set table to find the intersecting edges of T i with S i 12: if There are exactly 2 intersecting edges then 13: Locate the 2 intersections via univariate interpolation with respect to x 14: Sort the 2 resulting points by their polar-coordinate θ 15: Store the line-segment into R L 16: end if 17: Query the minimum distance to the AABB tree object, store the result into D 7: end for 8: Perform parallel reduction to find the minimum distance from D 9: end for 10: return all computed minimum distances.
Procedure 2: Parallel distance querying. 5 International Journal of Aerospace Engineering (4) Perform level set evolution, compute the gas generation rate to be used during the next iteration 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 expresses B 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 T E x of this probe point. Considering that the main axis of the SRM is aligned with one of the grid axes, the T E values fetched from a slice of the grid are equivalents of the T E values fetched from a sampling section. Like uniform models, the minimum one among the T E values fetched from a slice is taken as the final output of the nonuniform burnback model.

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: where L d S i is the thickness function of the designs to be judged, L e S i is the actual minimum required thickness to safely resist heat, and L d ′ S i 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: Input: D x , B for each sample point, p t data of the SRM 1: Discretize p t into a sequence of length n p such that ∀i < n p , P i = p t i 2: Build the propagation distance data sequence D p of length n p 3: ∀i < n p , assign D p|1 = t 1 r P 1 , D p|i = D p|i−1 + t i r P i 4: for all sampling sections S i parallel do 5: In sequence D p , find the first i such that D p|i−1 < D S i , B , D p|i ≥ D S i , B 6: Perform a linear interpolation within interval t i−1 , t i , find T E S i such that D p T E S i = D S i , B , 7: end for 8: return For all the sample points, T end − T E S i is the exposure duration. 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 Section 4.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.  Figure 4, 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.
The 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 2015 mm/MPa n . 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).
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: where T l is the duration before the carbonaceous char layer starts to move inward. The constants measured from the experiment data are T l = 1 3 s and R c R = 0 35 mm/s. 6.50e + 06 6.53e + 06 6.57e + 06 6.60e + 06 6.63e + 06 6.67e + 06 6.70e + 06 6.73e + 06 6.76e + 06 6.80e + 06 6.83e + 06 6.86e + 06 6.90e + 06 6.94e + 06 Other erosion models can be simply integrated by replacing the R c R T end − T E x item with corresponding numerical integrations.

Burnback Model.
The SRM is then fed to the above-implemented code. In Figure 7, 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.
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 Figure 8, 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.

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 L 0 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 L 0 ) 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).
From the tail dome of an actual SRM sample which had experienced a thrust experiment, we have measured the   International Journal of Aerospace Engineering actual eroded thickness of the HIL from 3 points (see Figure 10, 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. Using 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 − 4 m 3 and 1 35e − 4 m 3 , respectively. The design produced by the proposed approach is 3.85% lighter compared to that of the traditional approach.

Effect of Discretizing
Resolutions. The sampling resolution used in the above process is n s = 1024, n r = 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 n s and n r combinations and various SRM models. The results can be summarized by the following conclusions: (1) If H is a revolving surface, n r ≥ 5N as usually provides enough accuracy, where N as is the number of axial sweeping geometric features (e.g., fins); setting a larger n r while leaving n s unmoved produces little improvement in accuracy (2) If H is not a revolving surface, n r should be set such that the spacing between neighboring sampling points is smaller than the minimum characteristic scale of H A C B Figure 10: Positions of the 3 measured points.    Figure 12: A sample grain to verify the nonuniform burnback model. 9 International Journal of Aerospace Engineering (3) n s 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 Figure 12, 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. According to [21] and the ellipse equation, the following analytical result can be derived (assuming l ≫ a − h): where r s and r f are the constant burning rates of the slowerand faster-burning propellants, respectively.
Letting l = 100 mm, a = 40 mm, b = 15 mm, h = 15 mm, r s = 6 mm/s, and r f = 9 mm/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. Profiling 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.  Figure 4 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 n r = 180, n s = 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.
Obviously, the efficiency of the proposed algorithm is significantly higher than that of the AABB tree approach, with whatever discretizing resolution. More precisely, Procedure 1 saves roughly 35% of the execution time compared to the AABB tree approach.

Parallel Generating of Sampling Points: Procedure 1.
Here, we measure the improvement introduced by the parallel computing in Procedure 1. The triangular model having 4635 triangles is employed as the testing model. Accuracy parameters are set as n r = 180 and n s = 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.
In Figure 15, the numbered legends denote the n s 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 Distance querying rate (Mspps) Figure 16: Parallel distance querying profiling data. 11 International Journal of Aerospace Engineering curve slightly shifts upward as n s 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: Procedure 2. Using the same geometric model and hardware platform as Section 4.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 n s 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.

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

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

Conflicts of Interest
The authors declare that they have no conflicts of interest.