International Journal of Aerospace Engineering
Volume 2019 (2019), Article ID 5789430, 13 pages
https://doi.org/10.1155/2019/5789430

Research Article

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

Ran Wei, Futing Bao, Yang Liu, and Weihua Hui

Science and Technology on Combustion, Internal Flow and Thermo-structure Laboratory, School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China

Correspondence should be addressed to Weihua Hui; zhongyuancao@163.com

Received 27 July 2018; Revised 25 October 2018; Accepted 12 November 2018; Published 21 January 2019

Academic Editor: Marco Pizzarelli

Copyright © 2019 Ran Wei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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

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

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 informationP.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 unnecessaryP.3.There is no hole or flame-retardant material inside the grainP.4.The thickness of the HIL on the same radial section is uniform. Otherwise, manufacturing of the HIL and the grain would become difficultP.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 , and the set of the triangles on the inner surface of HIL as .

Letting be the duration before the flame arriving point , we have

where is the shortest flame traveling path from B to and is the transient burning rate at the burning location on . Since the burning rate has spatial dependency, the shortest path 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, is always the shortest line segment connecting and , so the solution of equation (1) can be greatly simplified (see the discussion in Section 2.2).

The burning rate is decided by the widely known exponential burning rate formula [15]:

where and are decided by the local propellant material property and initial temperature and is the current gas pressure near . Usually, the values of and can be determined by burning rate experiments and are constant among different SRMS, while 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 be the elapsed time after exposure, the erosion model at can be summarized by the following time-dependent function:

where denotes the recession rate at time , and denote the heat transfer caused by conduction and advection, respectively, denotes the mechanical effect of condensed particles, and stands for the material property of the HIL. The last parameter implies that function has a time-memory property, i.e., may output different values at different moments even if , , and remain identical.

Let be the moment when the SRM stops working. The thickness of the eroded HIL can be computed as follows:

In order to prevent the overheating of SRM cases, the thickness of the residual HIL at must not be too small. The ideal residual thickness can be decided according to the material properties of the SRM. Letting 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 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 is the length of the shortest line segment connecting to . 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 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 which has to be computed(2)Compute values of by the uniform or nonuniform burnback model(3)For each of the points, compute the required HIL thickness by the chosen erosion model

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 is enough to guide the design. We thus discretize the HIL into radial sections, and take sample points from each of the sections. and 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 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.

The above section-by-section discretizing scheme allows us to further simplify the solution process. For each section , the ideal HIL thickness equals the maximum value obtained from all its belonging sample points. Letting the set of all sampling points taken from be , the minimum allowed thickness on can be calculated as follows:

Since is a nonnegative function, monotonically decreases as grows. Consequently, equation (7) can be further simplified:

Compared to equation (7), equation (8) computes 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

Letting be the coordinate range of , the equation of the th radial section is , where . The intersection of and 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 is an infinite plane perpendicular to the axis. The and items can therefore be eliminated from the plane equation. In this situation, the intersection of any line and can be located via three independent univariate interpolations with respect to , thus evading the 3 × 3 matrix operation in general approaches(2)The facets in generally distribute perpendicular to , so a large portion of the facets can be excluded from the intersecting computation according to their 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 to can be safely ignored. In general methods, there have to be some extra operations to handle this case

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 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 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 of for each sample point are obtained. In order to solve equation (6), the interior ballistic data of the SRM, i.e., 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 and 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 iteration

Input: Triangle set , axial discretizing resolution , circumferential discretizing resolution

1: , compute its coordinate range

2: Sort all triangles in according to

3: Initialize resulting sampling point set

4:

5: for all plane at coordinate serial do

6:    Initialize resulting line-segment set

7:    Remove from H

8:    for all parallel do

9:        if then

10:            For each vertex of , compute its location identifier

11:            Look up in a pre-set table to find the intersecting edges of with

12:            if There are exactly 2 intersecting edges then

13:                Locate the 2 intersections via univariate interpolation with respect to

14:                Sort the 2 resulting points by their polar-coordinate

15:                Store the line-segment into

16:            end if

17:        end if

18:    end for

19:    Sort all line-segments in according to the polar-coordinate

20:    Compute the polar-coordinate range of all line-segments in

21:    On , generate rays which start from the center and distribute within 

22:    for all parallel do

23:        Fetch intersecting rays of according to

24:        Locate the sampling points by solving the 2D intersection equations

25:        Store the sampling points into

26:    end for

27: end for

28: return .

Procedure 1: Parallel sampling point generation.

Input: Triangle set , sampling point set

1: Build a AABB tree object from

2: for all sampling sections serial do

3:    Fetch all sampling points on from

4:    Initialize resulting array

5:    for all fetched sampling points parallel do

6:        Query the minimum distance to the AABB tree object, store the result into

7:    end for

8:    Perform parallel reduction to find the minimum distance from

9: end for

10: return all computed minimum distances.

Procedure 2: Parallel distance querying.

Input: for each sample point, data of the SRM

1: Discretize into a sequence of length such that

2: Build the propagation distance data sequence of length

3: , assign ,

4: for all sampling sections parallel do

5:    In sequence , find the first such that ,

6:    Perform a linear interpolation within interval , find such that ,

7: end for

8: return For all the sample points, is the exposure duration.

Procedure 3: Compute 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 expresses and 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 of this probe point. Considering that the main axis of the SRM is aligned with one of the grid axes, the values fetched from a slice of the grid are equivalents of the values fetched from a sampling section. Like uniform models, the minimum one among the 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:

where is the thickness function of the designs to be judged, is the actual minimum required thickness to safely resist heat, and 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:

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.

4.1. Uniform Approach

4.1.1. Sample SRM

The sample SRM is demonstrated in 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.

(a) The sectional view

(b) The initial burn surface

Figure 4: The grain of the sample SRM.

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 , . 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:

where is the duration before the carbonaceous char layer starts to move inward. The constants measured from the experiment data are and .

Substituting equation (11) into equation (4) and equation (5), we have the following compact equation to design the thickness of the HIL:

Other erosion models can be simply integrated by replacing the item with corresponding numerical integrations.

4.1.2. Burnback Model

The SRM is then fed to the above-implemented code. In Figure 7, the generated sampling points, colored according to their values, are demonstrated. The value and exposure duration on each sampling section is plotted in Figure 8 by the blue line and the red line, respectively.

Figure 7: on a sampling point set.

Figure 8: and exposure duration on sampling sections ( and ).

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 values from the CAD platform. It can be seen that every 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 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 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 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 ) 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: Resulting 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 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.

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

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 and , 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 is , , which produces an accurate 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 and combinations and various SRM models. The results can be summarized by the following conclusions:

(1)If is a revolving surface, usually provides enough accuracy, where is the number of axial sweeping geometric features (e.g., fins); setting a larger while leaving unmoved produces little improvement in accuracy(2)If is not a revolving surface, should be set such that the spacing between neighboring sampling points is smaller than the minimum characteristic scale of (3) 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 and , 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 ):

where and are the constant burning rates of the slower- and faster-burning propellants, respectively.

Letting , , , , , and , 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: 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 (Section 3.2.1)(2)The parallel computing technology utilized in Procedure 1(3)The parallel computing technology utilized in Procedure 2(4)The parallel computing technology utilized in the nonuniform burnback model

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.

4.3.1. Efficient Sampling Point-Generating Process: Procedure 1

Discretizing the CAD model shown by 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 , , 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, Procedure 1 saves roughly 35% of the execution time compared to the AABB tree approach.

4.3.2. 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 and . 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 Figure 15, the numbered legends denote the 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 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 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.

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.

References

Y. H. Xu, X. Hu, Y. X. Yang, Z. X. Zeng, and C. B. Hu, “Dynamic simulation of insulation material ablation process in solid propellant rocket motor,” Journal of Aerospace Engineering, vol. 28, no. 5, article 04014118, 2015.

M. Storti, L. A. Crivelli, and S. R. Idelsohn, “An efficient tangent scheme for solving phase-change problems,” Computer Methods in Applied Mechanics and Engineering, vol. 66, no. 1, pp. 65–86, 1988.

S. S. Katte, S. K. Das, and S. P. Venkateshan, “Two-dimensional ablation in cylindrical geometry,” Journal of Thermophysics and Heat Transfer, vol. 14, no. 4, pp. 548–556, 2000.

A. L. Murray and G. W. Russell, “Coupled aeroheating/ablation analysis for missile configurations,” Journal of Spacecraft and Rockets, vol. 39, no. 4, pp. 501–508, 2002.

Y.-K. Chen and F. Milos, “Three-dimensional ablation and thermal response simulation system,” in 38th AIAA Thermophysics Conference, American Institute of Aeronautics and Astronautics, Toronto, Ontario, Canada, June 2005.

B.-C. Yang, F.-B. Cheung, and J. Koo, “Prediction of thermomechanical erosion of high-temperature ablatives in the SSRM facility,” in 33rd Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, Reno, NV, USA, January 1995.

H. Wirzberger and S. Yaniv, “Prediction of erosion in a solid rocket motor by alumina particles,” in 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, American Institute of Aeronautics and Astronautics, Tucson, AZ, USA, July 2005.

L. I. Li, T. Yang, X. H. Cheng, and L. Yang, “Ablation model of silicon insulator in ramjet combustion chamber,” Journal of Propulsion Technology, vol. 33, no. 3, pp. 451–454, 2012.

M. A. Willcox, M. Q. Brewster, K. C. Tang, and D. S. Stewart, “Solid propellant grain design and burnback simulation using a minimum distance function,” Journal of Propulsion and Power, vol. 23, no. 2, pp. 465–475, 2007.

L. Yang, B. Fu-ting, C. Qiang, H. Haifeng, and D. Jiajia, “GPU computing architecture-based discrete voxelized burning surface calculation method for complex SRM grains,” Journal of Solid Rocket Technology, vol. 34, no. 1, pp. 18–22, 2011.

W. Sullwald, G. J. Smit, A. J. Steenkamp, and C. W. Rousseau, “Solid rocket motor grain burn back analysis using level set methods and Monte Carlo volume integration,” in 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, American Institute of Aeronautics and Astronautics (AIAA), San Jose, CA, USA, July 2013.

R. Wei, F. Bao, Y. Liu, and W. Hui, “Combined acceleration methods for solid rocket motor grain burnback simulation based on the level set method,” International Journal of Aerospace Engineering, vol. 2018, 12 pages, 2018.

K. Albarado, A. Shelton, and R. Hartfield, “SRM simulation using the level set method and higher order integration schemes,” in 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, American Institute of Aeronautics and Astronautics (AIAA), Atlanta, GA, USA, July 2012.

K. Toker, H. Tinaztepe, and M. Aksel, “Three-dimensional internal ballistic analysis by fast marching method applied to propellant grain burnback,” in 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, American Institute of Aeronautics and Astronautics, Tucson, AZ, USA, July 2005.

A. I. Atwood, T. L. Boggs, P. O. Curran et al., “Burning rate of solid propellant ingredients, part 1: pressure and initial temperature effects,” Journal of Propulsion and Power, vol. 15, no. 6, pp. 740–747, 1999.

K. E. Hoff, J. Keyser, M. Lin, D. Manocha, and T. Culver, “Fast computation of generalized Voronoi diagrams using graphics hardware,” in Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques—SIGGRAPH ‘99, pp. 277–286, Association for Computing Machinery (ACM), New York, NY, USA, 1999.

A. Sud, M. A. Otaduy, and D. Manocha, “DiFi: fast 3D distance field computation using graphics hardware,” Computer Graphics Forum, vol. 23, no. 3, pp. 557–566, 2004.

P. Alliez, S. Tayeb, and C. Wormser, “3D fast intersection and distance computation,” in CGAL User and Reference Manual, CGAL Editorial Board, 4.13 edition, 2018.

M. de Berg, Computational Geometry: Algorithms and Applications, Springer, Berlin; New York, 2nd edition, 2000.

J. Reinders, Intel Threading Building Blocks, O’Reilly & Associates, Inc, Sebastopol, CA, USA, 1st edition, 2007.

X. Lin, B. Wang, and S. Jin, “A method for the interior ballistics calculations of solid rocket motors with dual burning rate propellant,” Journal of Solid Rocket Technology, vol. 4, pp. 4–4, 1991.