Solid rocket motor internal ballistics using an enhanced surface-vorticity panel technique

It may be safely argued that the tools that are available for characterizing the internal gas dynamics in the context of solid propellant rocket motor design have evolved into three fundamental groups. These consist of (i) low-order and mainly one-dimensional approaches that return bulk flow values; (ii) analytical solutions for special cases; and (iii) high-resolution computational fluid dynamics (CFD) solutions. On the one hand, low-order solutions tend to be of limited fidelity and utility for design considerations entailing performance enhancement, regression-rate control, acoustic instability mitigation, and so on. On the other hand, high-resolution CFD solutions remain of limited value in design iteration activities due to their steep computational cost. Clearly, the need exists for a formulation that can provide detailed solutions with tractable efficiency. The main purpose of this work is to present indeed such a method by leveraging the characteristics and versatility of an enhanced surface-vorticity (SV) framework.

The more detailed solutions of the internal ballistics problem for solid rocket motors offer a foundation upon which erosive burning models, performance prediction paradigms, and acoustic stability platforms for rocket motor analysis and characterization can be undertaken. The simplest departure from bulk mode considerations is a one-dimensional solution outlined with lucidity by Barrere et al.¹ The basic assumptions for this solution are comparable to those underlying potential flow theory. In fact, it has been demonstrated by several researchers that similar potential flow assumptions can be applied to develop analytical solutions to rocket motor internal flows with remarkable fidelity. Examples include those by McClure et al.,² Hart and McClure,^3,4 Yagodkin,⁵ Terrill and Colgan,⁶ Haloulakos,⁷ Smith-Kent et al.,⁸ Salita,⁹ Majdalani and Akiki,¹⁰ Kurdyumov,¹¹ Saad and Majdalani,¹² Maicke et al.,¹³ Majdalani and Saad,¹⁴ and Williams and Majdalani¹⁵ to name a few.

As affirmed by Culick,¹⁶ modeling internal flowfields proves to be essential in propulsive applications, especially those that are based on solid rocket motors. The instantaneous flowfield has a strong bearing on prescribing the thrust and ballistic performance of the motor at hand, as it directly affects the behavior of pressure oscillations, grain regression, particle–mean flow interactions, erosive burning, nozzle erosion, slag formation, and flight rotation.

For example, in several dedicated studies by Apte and Yang,^17–20 Fabignon et al.,²¹ Flandro et al.,²² Chedevergne et al.,²³ Boyer et al.,^24,25 and others, the instantaneous field is decomposed into a steady bulk motion and one or several unsteady disturbance modes. In practice, the bulk motion refers to the average, steady-state gaseous flow structure, which is routinely represented by the solution of a wall-bounded porous duct (channel or tube) with distributed wall injection. As for the oscillatory field, it is usually associated with unsteady disturbances that are pressure-driven (curl-free, acoustic); boundary-driven (divergence-free, vortical); and entropy-generated (biharmonic, thermal). The reader may consult in this regard with Chu and Kovásznay²⁶ on mode splitting and Bhatia et al.²⁷ on the Helmholtz–Hodge decomposition principle. An excellent overview of decomposition techniques is also provided by Ewert and Schröder.²⁸ The time-dependent disturbances can also comprise hydrodynamic instability waves, which are triggered by the breakdown modes of the mean flow; these generally evolve at low frequencies than the vortico-acoustic and entropy modes and, as such, can be resolved independently (Varapaev and Yagodkin,²⁹ Beddini and Roberts,^30,31 Griffond et al.,³² and Gazanion et al.³³). They can also be resolved in combination with the entire set of unsteady disturbance mode structures, as shown by Majdalani et al.,³⁴ using a compressible hydrodynamic stability solver. More importantly, the necessity of rapidly predicting the mean flow structure for a given geometric configuration can be of profound interest given the strong coupling that exists between the steady and unsteady flowfields.

Although the earliest stability investigations of solid rocket motors can be seen to reduce the motor chamber to a porous enclosure, they are not always based on realistic conditions. In fact, the first theoretical mean flow model used in conjunction with rocket instability may be traced back to Grad.³⁵ Therein, Grad neglects the effect of mean flow by assuming a quiescent, stagnation flow medium inside a simulated rocket chamber; the resulting characterization remains limited to the so-called “organ-pipe” or Helmholtz-type oscillations.

By way of improvement, the work of Hart, McClure, and coworkers may have been perhaps the first to adopt a non-vanishing mean flow profile in the context of rocket flow analysis. Their work is carried out in an effort to resolve the thin region above the propellant surface where the most notable pressure and velocity coupling mechanisms take place (Hart and McClure;^3,4 Hart et al.;³⁶ Hart and Cantrell;³⁷ and McClure et al.³⁸). In actuality, the acoustic stability treatment of T-burners by McClure et al.² may have been the first to incorporate a non-vanishing bulk flow model. In this context, their solution corresponds to the ideal motion of a frictionless fluid sweeping down the bore of a porous cylinder. Despite its irrotational character, the so-called Hart–McClure profile is later shown to generate favorable predictions of motor performance, namely, by helping to identify the intrinsic coupling that occurs between the steady and time-dependent fluctuations.

Shortly thereafter, Culick³⁹ manages to derive another inviscid representation of the mean flow in a circular-port motor using a rotational framework. This solution is later verified experimentally by Dunlap et al.⁴⁰ At first glance, Culick's profile offers the advantage of satisfying the no-slip requirement at the sidewall, while allowing slip at the headwall. The profile itself proves to be identical to a formulation devised a decade earlier by Taylor,⁴¹ albeit in the context of a paper manufacturing application. For this reason, the ensuing model is often referred to as the “Taylor–Culick” profile, named as such after both scientists who arrived at the same formulation independently. In fact, a plethora of propulsion-related studies would later make extensive use of the Taylor–Culick model as the basis for their investigations. On this count, one may cite Beddini,⁴² Majdalani and Van Moorhem,⁴³ Majdalani,⁴⁴ Griffond et al.,³² Majdalani et al.,⁴⁵ Flandro and Majdalani,⁴⁶ Chedevergne et al.,^23,47 Gazanion et al.,³³ and Bouyges et al.⁴⁸ However, except for its ability to capture wall vorticity very well, many characteristics associated with the Taylor–Culick profile are adequately captured by the irrotational Hart–McClure model, especially when area averages are taken along the motor axis, as shown by Saad and Majdalani.¹² 

In extending the internal flow models beyond the Hart–McClure and Taylor–Culick profiles, several formulations have been advanced for simple geometric configurations. Some of these seek to compensate for the inability of the foregoing models to account for certain aspects of the flow, particularly, those that can become significant under specific conditions such as high-speed motion or wall regression. In this vein, one may explore the benefits of incorporating the effects of compressibility, wall regression, grain taper, viscosity, headwall injection, or swirl.

First, the effects of compressibility are addressed in several studies that include integral formulations by Traineau et al.,⁴⁹ Balakrishnan et al.,⁵⁰ Gany and Aharon,⁵¹ and Akiki and Majdalani.^52,53 These are nicely complemented by closed-form Rayleigh–Janzen expansions in the wall Mach number squared by Majdalani,⁵⁴ Maicke and Majdalani,⁵⁵ and Cecil and Majdalani.⁵⁶ 

Second, the sensitivity of the internal flowfield to wall regression or taper is examined assuming a planar configuration by Dauenhauer and Majdalani;^57,58 Majdalani and Zhou;⁵⁹ Boutros et al.;⁶⁰ Asghar et al.;⁶¹ Xu et al.;⁶² and, more recently, by Majdalani and Xuan.⁶³ In the latter, a viscous boundary layer correction is provided to overcome the singular pressure distribution and its normal gradients near the midsection plane of an expanding porous channel. As for the axisymmetric porous tube with wall regression, it is first treated by Uchida and Aoki,⁶⁴ and later re-examined by Majdalani et al.,^65,66 Boutros et al.,⁶⁷ and Saad and Majdalani.⁶⁸ It may be instructive to note that, in regression-sensitive studies, viscous effects near the wall are typically retained.

Third, in an effort to account for irregular grain perforations, several enhancements to the Hart–McClure and Taylor–Culick profiles are pursued, often with a significant penalty in the complexity of the resulting formulations. On this note, approximations have been obtained by Clayton,⁶⁹ Saad et al.,⁷⁰ and Sams et al.,⁷¹ who modified the internal-burning rocket chamber by allowing for small wall taper; other, more significant efforts include those by Kurdyumov,⁷² Van Horn and Majdalani,⁷³ and Bouyges et al.,⁴⁸ who extended the Taylor–Culick solution to chambers with irregular cross sections, such as those with star-shaped perforations.

Fourth, in moving closer to the central topic of this article, formulations are considered wherein the internal flowfield is reconstructed by allowing for arbitrary headwall injection.⁷⁴ The corresponding problem is analyzed in both axisymmetric and planar configurations by Majdalani and Saad,⁷⁵ and Saad and Majdalani,⁷⁶ respectively. In subsequent work, Saad and Majdalani¹² manage to identify solutions of the Taylor–Culick type with varying kinetic energies. Using a variational procedure anchored around Lagrangian multipliers, these researchers are able to unravel a wide range of motions that extend from strictly irrotational to highly rotational fields. Finally, the effects of a swirling motor case are addressed by Bastress,⁷⁷ Manda,⁷⁸ King,⁷⁹ Lucy,⁸⁰ Dunlap,⁸¹ Johnson and Lecuyer,⁸² Fist and Majdalani,⁸³ and Cecil and Majdalani.⁸⁴ Work in the latter leads to two distinct analytical formulations that match the finite volume numerical simulations of a spinning circular-port motor virtually identically.

Let us now recall that the assumptions of reversible and irrotational motion constitute the foundational requirements used to develop potential flow solutions for external motions. In the analysis of internal motions, such as those that apply to wall-bounded porous tubes, potential functions are often used; however, for external aerodynamic analysis, the doublet or vortex filament approach is more frequently used because of the extension of the distributed circulation to the Kutta–Joukowski theorem; the latter relies on circulation in the calculation of lift and induced drag. Several different implementations of this approach to aerodynamics problems have, in fact, led to the creation of so-called “panel codes.” These include, but are not limited to, Pan Air; PMARC; CBAERO; Vorlax; VSAero; VSPAero; and, most recently, FlightStream^®.⁸⁵ The latter, whose framework constitutes the underlying vehicle of this study, brings the speed and fidelity of potential methods to bear on the solid rocket motor internal ballistics problem; this is accomplished while offering an intuitive user-friendly interface, flow separation modeling, and a high degree of interactivity with other modern engineering tools such as computer aided design (CAD) and finite-element stress analysis packages. In this work, it is shown that, through reliance on a surface-mesh-based numerical approach to extend the foundational potential theory used to develop analytical solutions to solid rocket motors, high speed computations with fidelity suitable for physical interpretation and engineering design can be expanded to a much broader range of geometric configurations. Through this approach, many of the complexities facing further developments of analytical models for rocket motors can be overcome, thus providing a comparably straightforward and accurate formulation of the problem in the case of irregular grain geometry.

In this article, the computational model is introduced using a description of the geometric creation process, the meshing process, and the flow solver characteristics. This is followed by several verification cases that explore the ability of the solver to produce accurate solutions to basic theoretical problems for which analytical formulations can be found in the literature. These include the fundamental cases of an internal-burning cylinder and a slab motor with and without wall taper and headwall injection. A validation study is also pursued by comparing the simulation results of the present framework to full-scale computational studies as well as experimental measurements acquired for the reusable solid rocket motor (RSRM), thus illustrating the degree of fidelity as well as the underlying limitations of a surface-vorticity solver.

All of the models in this study are originally constructed in a computer-aided design (CAD) software (e.g., SolidWorks) and imported into the surface-vorticity FlightStream^® solver as initial graphics exchange specification (IGES) surfaces. The surfaces correspond to the exposed boundaries of the solid propellant grain. In fact, the solver features an intuitive user-interface and a CAD module with several practical tools to assist in generating a proper surface mesh. It also contains an assortment of built-in meshing tools to improve raw mesh imports including trimmed and aligned meshers, which have been used extensively throughout this study. The trimmed mesher, which is generally applied to irregularly shaped surfaces, generates unstructured isotropic facets. The aligned mesher creates structured anisotropic facets, with the latter being more effective on elongated surfaces that require a more refined mesh in one particular dimension vs another; these anisotropic facets become quite convenient to implement in modeling the long tubular sections of a solid rocket motor.

It proves helpful to ensure, at the outset, that the mesh surface appears as a single continuous structure because the surface-mesh solvers are not overset. The solver exhibits a built-in diagnostic tool as well as a variety of repair tools to correct these meshing errors, should they occur. It is also desirable to ensure that the mesh quality remains consistently favorable with respect to the facet aspect ratio. In extreme cases, it may be beneficial to simplify the modeled geometry to avoid the pitfalls of poor quality meshes.

Within the solver, the burning surfaces can be declared to be Dirichlet inlet boundaries that impose a nonzero velocity on the surface of the mesh by taking into account the propellant burn rate. The surface normal specification of the inlet boundary is widely adopted in modeling solid propellant grains because it applies a realistic injection or blowing velocity along the surface in the normal direction of each mesh facet. On this note, the model requires setting the surface normal unit vectors of all mesh faces facing inward. To assist with this operation, an option can be toggled in the solver to colorize the faces based on the direction of their surface normals relative to their current view. In addition, the surface normals can be flipped rather easily with a built-in toggle. From a practical perspective, the surface-vorticity solver enables the user to specify a custom inlet profile should a variable inlet velocity be desired. Naturally, this feature becomes valuable for modeling erosive burn rate augmentation and non-homogeneous propellant formulations, among other effects. This option accepts as input a text file of coordinates and corresponding inlet velocities. The specified inlet velocities can then be applied to all mesh surfaces that are adjacent to the coordinates prescribed in the text file. In most cases, it is helpful to specify a custom profile down the centerline of the chamber cavity. This custom inlet feature has been developed and integrated into the surface-vorticity solver to enhance its capabilities and effectiveness at resolving internal flows such as those encountered in rocket motors.

Our principal simulation tool consists of a high-order surface-vorticity solver that assigns vortex rings to each mesh facet, where each edge is treated as a segmented component of the vortex. These vortex rings interconnect, sharing edges with neighboring rings, and encompassing the entire surface body of the input geometry. The induced velocity is calculated for each segment of these vortex rings according to the Biot–Savart law, specifically depending on the plane created by the segment's vertices and the spatial point of induction. The induced velocity of a particular vortex ring is found through summation of the induced velocities of each segment forming the ring. In this process, the direction of each vortex ring is determined by the facet's normal vector and, as such, it is imperative to configure the facet normals homogenously. The vortex strength is acquired from known velocities at the boundaries in accordance with the Helmholtz law. To overcome the inversion errors that can occur during the manipulation of the large matrices that accompany large surface meshes, the mesh may be subdivided into parts that can be individually solved in an iterative manner. Consequently, the solver's convergence rate becomes dependent on the stability of the surface vorticity and wake integration plane.⁸⁵ 

Inherent to the vortex formulation of the potential flow equations, a freestream velocity must be defined. In extending our mathematical framework to internal flows, this freestream velocity becomes immaterial and can be set to an arbitrarily small value. As for the definition of the solid rocket motor problem, it begins at the fluid flow level. Instead of representing a regressing surface, mesh surfaces can be used to capture the inception of gaseous motion. In practice, the gaseous injection velocity can be readily calculated using mass conservation that takes into account the propellant burn rate and the corresponding conversion of the solid propellant into a gaseous substance.

Finally, it may be instructive to note that the surface vortices in the outflow section of a rocket chamber can become unbounded, thus rendering the solution unphysical in the exit plane. In the present work, this deficiency is readily overcome by extending the length of the chamber in a manner to ensure that the unphysical outlet falls well beyond the domain of investigation. In all of the test cases considered here, a convergence threshold of 10⁻⁶, which is representative of the stability of the surface vorticity, is set as a requirement. This threshold is found to be fairly adequate in producing accurate internal flow solutions. As a result, smooth and direct solution convergence will be reported for all of the tested models.

In order to provide a basis for theoretical comparisons to the surface-vorticity solutions, the analysis starts by deriving irrotational mean flow profiles for the basic cases of cylindrical and planar chambers with slight taper and uniform headwall injection; this effort may be guided by the works of Saad et al.⁷⁰ and Sams et al.,⁷¹ where rotational equivalents are developed. In the interest of clarity, a brief overview of the resulting formulations is presented here; for a more in-depth description of the underlying procedures, the aforementioned studies may be referred to.

Figure 1 illustrates the control volume and principal parameters that may be used to model the axisymmetric, circular-port, cylindrical case with a small taper angle. The axis of symmetry is anchored at the center of the headwall while the burning surface is oriented at a slight taper angle $α$⁠. The surface injection speed is defined to be $ub$⁠, where the subscript “b” refers to either a “burning surface” or a “blowing” property. Along similar lines, a uniform headwall injection speed is taken to be $uh$⁠, where the subscript “h” reflects a headwall property.

To make further headway, it is convenient to non-dimensionalize all variables of interest. For this purpose, the burning surface injection velocity $ub$ and the initial radius $r0$ of the propellant burning radius $rb$ are chosen as reference values. Then, using overbars to denote dimensional quantities, one can set

where $Ψ$ alludes to the mean flow normalized streamfunction; the latter is prescribed by the incompressible Stokes equation, which comprises a function $f(ψ)$ that depends on the vorticity at the boundaries; one gets:³⁹ 

For a potential motion, one sets $f(ψ)=0$ such that Eq. (3) simplifies into a linear partial differential equation for steady, incompressible, irrotational, and axisymmetric flow viz.

The normalized boundary conditions become

These physical conditions ensure that no flow crosses the centerline; they also specify the ratio $β$ of the headwall injection speed to the burning surface injection speed as well as the two normalized components of the injection velocity at the burning surface. Given that the radial velocity in a non-tapered chamber remains strictly one-dimensional, it is reasonable to assume similar behavior for slightly tapered chambers. This particular assumption enables us to simplify the corresponding Stokes-like equation^70,71 to an ordinary differential equation that can be solved rather straightforwardly. Assuming steady, irrotational, axisymmetric, and incompressible flow with a small taper angle, one can set

This permits reducing Eq. (4) into:

To solve the reduced equation of motion, congruent boundary conditions must be devised for the streamfunction; one may follow Sams et al.⁷¹ and write

A suitable streamfunction can be determined using a twofold integration to obtain

At this stage, the velocity components can be readily determined through partial differentiation of Eq. (9). Using Stokes's streamfunction definition, one gets

and

The pressure follows directly from Euler's equation. Given that the solution is irrotational, the latter may be integrated between any two points, thus leading to Bernoulli's expression. By selecting the burning surface b as a reference point, and recalling that $ρg$ represents the gaseous density, one can put

and so

Then using $qb=12ρgub2$ for the surface dynamic pressure, the pressure coefficient can be extracted directly from

The problem can be reframed for the case of a planar slab motor chamber with porous walls where, instead of axisymmetric flow, no crossflow may be imposed. Figure 1 can just as fittingly represent the planar geometry, using $x$ for the horizontal abscissa and $y$ for the vertical ordinate. The normalized boundary conditions for this configuration translate into the following:

As before, the small angle approximation can be made to simplify the potential flow equation in Cartesian coordinates. Assuming steady, irrotational, and incompressible flow with a small taper angle and no crossflow, the potential equation of motion collapses into

Following a similar procedure to that used by Saad et al.,⁷⁰ a set of two boundary conditions may be specified, namely:

These could be implemented in conjunction with a two-step integration of Eq. (16) to retrieve

Here too, the velocity components may be readily determined from the planar streamfunction, specifically

and

As for the pressure coefficient, it may be extracted from Bernoulli's expression relative to the burning surface

At this juncture, practical limits that are prescribed by mass conservation⁷¹ can be imposed on the maximum headwall injection velocity. Since the flow speed converges to a constant value in a sufficiently long tapered channel,⁷¹ the analytical mean flow solution will remain physical so long as the headwall injection velocity does not exceed this asymptotic limit, $βlimit$⁠. For both axisymmetric and planar configurations, the headwall injection velocity ratio must be taken such that^70,71

Given that the models in question entertain modest taper angles, a small angle approximation is appropriate. In all cases considered, the headwall injection velocity ratio is found to fall well below $βlimit.$

To set the stage, mesh surfaces for both non-tapered and 2° tapered cylindrical chambers are shown in Fig. 2. Note that axial dimensions are provided to distinguish between the domain of investigation and the excess chamber length that is added in the downstream direction to promote a smooth outflow. The total area of each surface is declared as an inlet boundary. A standard structured grid is implemented on the headwall of the non-tapered chamber while a radially structured grid is applied to the headwall of the tapered chamber. In actuality, and after several simulation runs, no discernible differences are found between the solutions based on these two meshing stencils.

To illustrate the agreement between the surface-vorticity solver's predictions and the analytical solution given by Eq. (9), non-dimensional velocity contours are provided in Fig. 3 for the cylindrical configurations. In these two cases, the headwall injection speed is suppressed in order to better examine the effects of sidewall injection. It may be clearly seen that the differences between the surface-vorticity predictions (represented by full lines) and the analytical solution (broken lines) given by Eq. (9) remain virtually imperceptible except near the boundaries in both cases considered. One may also infer the essentially one-dimensional nature of the irrotational motion. This behavior is distinctly different from that of rotational flow solutions wherein the no-slip condition is secured along the burning surface and the mass flow remains much more concentrated near the central axis. In this comparison, the error increases near the mesh surfaces, which happens to be a characteristic of surface-vorticity methods: near-surface measurements are skewed by the singularities that accompany the vortex rings forming the mesh surface. In fact, work is presently underway toward the development of a vortex-splitting method to attenuate the corresponding error.

Besides this cursory comparison of velocity contours, attention may be turned to the streamline patterns. The latter can be generated for the two cylindrical models with various headwall injection velocities. To do so, the graphing tool's built-in streamline plotting feature is utilized;⁸⁶ this takes as input radial and axial velocity matrices and produces streamlines for a specified stream-field density. Forthwith, the analytical and numerical velocity fields are produced and plotted with equivalent densities in Fig. 4. Note that the surface-vorticity results are shown with full lines while the analytical model is depicted using broken lines. Since the graphing tool does not allow for a nonrectangular domain in the generation of streamline plots, the radial distance is normalized by the burning surface radius in the case of the tapered motor. For each of the geometric configurations, three headwall injection velocity ratios are examined: a zero headwall injection, a unitary injection ratio, and an injection ratio of 5. These are featured on the graph from top to bottom, respectively.

Unsurprisingly, the slip condition can be observed in the streamlines emerging from the burning surface, i.e., consistently with the irrotational flow assumptions. Further down the chamber, the tangential component of the flow at entry becomes even more appreciable. The streamlines lose curvature as the headwall injection speed is increased, thus causing the axial streamtube orientation to become more dominant. In both the tapered and non-tapered cases, the surface-vorticity streamlines computed for all three headwall injection ratios align quite convincingly with the analytical model.

To make further headway, one may take into account the absence of mean flow rotationality and proceed to deduce the pressure distribution directly from Bernoulli's energy conservation principle. Specifically, one may determine the change in centerline pressure over a range of taper angles using Eq. (14) between a reference point at the burning surface and any other point along the centerline. One gets

Note that the headwall injection velocity is relaxed in displaying these results. As shown in Fig. 5, where the pressure coefficients are compared for taper angles of $α=0°$⁠, $1°$⁠, $2°$⁠, and $3°$⁠, the surface-vorticity solution is found to be nearly identical to the analytical solution for all of the tested cases. These observations help to highlight the impact that even small taper angles can have on the flowfield. Although one-dimensional approaches typically overpredict the pressure drop in chambers with tapered aft ends, one can achieve quite reliable predictions using these simple models.

The same procedure may be repeated for the planar, slab motor case geometry. Using a similar gridding strategy to the one employed previously, the surface-vorticity mesh surfaces for the non-tapered and 2° tapered planar chambers are constructed and displayed in Fig. 6. For these models, it is important that the inert sidewalls are set at a sufficient width of at least six half heights to prevent edge effects in the spanwise direction from unnecessarily disturbing the two-dimensional flowfield.⁸⁷ 

The non-dimensional velocity contours and streamline patterns for the non-tapered and 2° tapered planar geometries with no headwall injection are presented in Figs. 7 and 8, respectively. With no headwall injection and a relatively small burning surface injection magnitude (as is typically the case in solid rocket motors), the strength of the flowfield is approximately half of that for a cylindrical chamber of similar dimensions. The surface-vorticity predictions and analytical formulation given by Eq. (18) show very good agreement for both cases in Fig. 7.

Here too, the error near the burning surface is expected due to the vortex singularities that are inherent to the surface-vorticity approach. The error appears slightly larger than in the cylindrical chamber, though it is likely just larger relative to the bulk flow magnitude. Moreover, as the axial distance increases, the surface-vorticity computations begin to lag behind the analytical solution, especially in the tapered case. The error becomes more concentrated near the midsection plane and grows in magnitude as the axial distance increases. At first glance, one may question whether this behavior could be attributed to the finite width of the models within the surface-vorticity solver. To test this hypothesis, wider models are constructed to determine if a larger spanwise dimension will help to alleviate this error. Forthwith, the solutions from the wider models are found to contain smaller errors; the observed error must be attributable to some other underlying issue that remains unidentified.

The streamline patterns in Fig. 8 are generated in the same manner as in the cylindrical case, i.e., using the streamline feature in the available graphing tool⁸⁶ for three headwall injection velocity ratios. The surface-vorticity results are presented in Fig. 8 using full lines whereas the analytical model is displayed using broken lines. Again, for the 2° tapered geometry, the radial distance is normalized by the burning surface radius $yb$⁠. At the outset, the surface-vorticity streamlines are seen to agree with the analytical model exceedingly well. Interestingly, the gradually increasing error seen in the contour plots does not appear to have an appreciable effect on the streamlines.

Finally, the midplane pressure profile is calculated for various taper angles using the same formulation of Bernoulli's principle given by Eq. (23). For these cases, the headwall injection velocity is excluded. The slight error observed in the velocity contour plots can also be seen here in Fig. 9, specifically increasing as the chamber is further elongated. The error appears to be somewhat uniform for all of the tested taper angles and is unlikely to be a function of the angle itself.

In seeking to achieve a more realistic and complex validation case, attention is turned to the space shuttle's reusable solid rocket motor (RSRM). The propellant's actual properties and geometric configuration may be retrieved from Thiokol's formal design proposal.⁸⁸ The model is based on motor operation past ignition transients, namely, when a fully developed flow is established assuming negligible changes to the initial geometry. In moving forward, a geometric CAD model of the exposed propellant surface is meticulously created in SolidWorks and then imported into the surface-vorticity solver as an IGES file; this surface is depicted in Fig. 10.

The RSRM consists of four segments: a forward finocyl segment; two central, slightly tapered (less than 1°) segments; and an aft segment that exhibits a larger taper angle at the exit. At the end of the tapered aft section, the chamber opens wider to allow the submerged nozzle to protrude into the chamber; this section as well as the nozzle are not modeled here. The finocyl consists of 11 fins with a taper at the outlet of the fins. Between each of the motor segments, radial slots may be identified and these become inhibited on the aft sides. Moreover, the first radial slot may be seen to be much smaller in size than the other slots. The entire surface is declared as an inlet boundary in the surface-vorticity solver except for the inhibited backsides of the radial slots, which are left as ordinary surfaces.

The mesh consists of approximately 28 500 facets; a close-up of the surface mesh is furnished in Fig. 11. A slight simplification of the geometry is necessary to improve the quality of the mesh at the outlet of the finocyl section. Therein, the fillets on the aft portion of the fins are removed to permit the use of the solver's aligned mesher on the central triangular faces. This change has a minor impact on the overall flowfield but is, nonetheless, necessary to ensure an accurate solution in this region.

The gaseous injectant is assumed to be calorically perfect, single phase, and ideal. In the interest of clarity, the propellant properties obtained from Thiokol's design proposal⁸⁸ are cataloged in Table I. Note that the mean chamber pressure is taken to be 6.21 MPa while the speed of sound $a0$ can be readily calculated based on the propellant flame temperature $Tf$ using

As for the solid-fuel regression rate $ṙ$⁠, it may be estimated using the empirical burn rate equation

where $a$ and $n$ stand for the temperature coefficient and pressure exponent, respectively. At this juncture, a relation between the regression rate and the injection velocity along the burning surface $ub$ can be deduced from mass conservation. Then, using $ρp$ and $ρg$ for the solid propellant and gaseous densities, one can write

These, along with other related flow properties, are calculated and documented in Table II. Before engaging in an RSRM validation test case, however, the potential flow solution produced by the surface-vorticity solver will be first verified through comparisons to one-dimensional flow and an isentropic process analysis of the streamtubes. The solution will then be compared to other computational models and openly available experimental data of RSRM qualification motors. Due to the high length-to-diameter and low port-to-throat ratio, the flow in the RSRM combustion chamber reaches relatively high Mach numbers. The flow is also known to be highly rotational due to the mass addition that is sustained along the burning surfaces as well as the effects of the radial slots between motor segments.⁸⁹ Thus, in order to judiciously carry out these comparisons, it is useful to implement basic rotational and compressibility corrections.

Our rotational correction is based on a relation between the irrotational analytical model of Sec. III A 1, given by Eq. (9), and the rotational equivalent derived by Sams et al.⁷¹ Considering that the RSRM is mostly composed of slightly tapered chambers, it is hoped that the analytical models will adequately represent the internal motion. As vorticity is generated at the wall during the injection process, the flow presumably develops until it reaches the maximum entropy state that is intrinsic to the rotational profile otherwise known as the Taylor–Culick profile.¹² For a tapered cylindrical chamber, the latter may be given by⁷¹ 

In both the irrotational and rotational models, the axial velocity proves to be linearly dependent on the axial distance; as such, a linear relation can be devised between them. At the centerline, the radial velocity of both models approaches zero. Along the centerline, the relations become

where the solution from the surface-vorticity solver, denoted by “SV,” represents the irrotational solution. Conversely, the subscript “ROT” is used to mark rotational solutions. From this exercise, one can readily infer that by multiplying the surface-vorticity's centerline velocity by $12π$⁠, the exact rotational velocity along the centerline may be recovered. A similar relation can be further confirmed from the planar analytical model, although it has no bearing on the RSRM analysis. This result is obtained owing to the attendant linearity with respect to the axial distance irrespective of whether headwall injection is present or not. In fact, Eq. (28) provides a relation between the maximum rotational velocity and the average bulk speed. The region that potentially deviates from this simple conversion factor is the central bore of the finocyl segment. Here, the flow will be three-dimensional, as the flow injects into the bore from the fin slots. This correction is therefore applied at the outlet of the finocyl and onward.

In rotational flow, one also recalls that the no-slip condition is upheld along the burning surface, thus resulting in a stagnation pressure loss near the surface. The gas propellant regression rate, gas density, and gas injection velocity are all functions of the pressure. As such, the injection velocity increases as the total pressure drops along the burning surface. In this analysis, one finds that the injection velocity is raised by nearly 10% toward the rear of the RSRM. To account for this, the increase in the burning surface injection rate may be iteratively coupled to the surface-vorticity solution through the use of a custom inlet profile as outlined next. First, an initial solution is found with constant injection speed throughout the enclosure. With this solution, the centerline static pressure may be directly calculated assuming isentropic flow along the centerline. The total pressure becomes essentially constant down the bore where the central streamline may be seen to originate at the headwall of the chamber. The pressure may then be calculated from

At any given axial station, the pressure of the injecting fluid at the burning surface is assumed to be approximately equal to the static pressure at the centerline. With this pressure profile, the new regression rates, injecting gas densities, and gas injection velocities may be straightforwardly determined. The recalculated gas injection velocity distribution may be fed back into the surface-vorticity solver, thus leading to a more accurate solution. In practice, the solution converges rather quickly, typically, after two or three iterations.

It is well known that relatively high Mach numbers are achieved in RSRM combustion chambers. To account for dilatational effects, the compressible Taylor–Culick profile derived by Majdalani⁵⁴ is presently used. The latter is derived using perturbation theory and is particularly useful because the first-order compressibility correction can be readily incorporated into the surface-vorticity solution in FlightStream^®, namely, after applying a rotational correction. This correction, denoted by $ζ$⁠, affects the axial velocity; the latter may be adjusted using

with

where

In the above, $Γ≈0.577 215 66$ represents Euler's gamma constant, while $κ0$ and $κ1$⁠, which appear in the compressible Taylor–Culick solution,⁵⁴ are evaluated and used in the surface-vorticity model. Note that the surface-vorticity solution, $(uz)SV,$ has been multiplied by $12π$⁠; this factor represents a one-to-one conversion multiplier of the irrotational solution to the rotational Taylor–Culick profile along the centerline. The correction is applied just aft of the finocyl section and is shifted to account for the flow velocity at the exit of the finocyl. The corresponding headwall velocity ratio used by Majdalani⁵⁴ is captured through $β$⁠. In the foregoing, the axial distance $z¯$ refers to the outlet of the finocyl, thus starting at a position of 4.5 m from the headwall. Finally, the average radius of the two central segments, $Ravg,$ may be employed as the normalizing radius in Majdalani's closed-form compressible correction.⁵⁴ 

It is important to note that this correction is derived for a non-tapered cylindrical chamber with a constant burning surface injection speed. As such, a small error may be expected to accompany this correction in the context of RSRM analysis. However, given that the correction appears at the first order in $Mw2$⁠, the corresponding error is neglected hereafter.

To begin, a comparison is carried out between a simple, continuity-satisfying, one-dimensional flow model and the surface-vorticity solution with no correction. In the foregoing analytical models, the axial velocity of the irrotational solution remains strictly invariant in the radial direction. Considering that the RSRM is mostly comprised of slightly tapered cylindrical chambers, a good agreement may be expected between the surface-vorticity solution and one-dimensional flow predictions. The one-dimensional solution may be calculated using conservation of mass at discrete axial stations throughout the motor. One simply gets

In the above, the inlet area $Ain$ represents the upstream burning surface and the outlet area $Aout$ corresponds to the cross-sectional area of the chamber at the specified axial location. The computed velocity is taken at the centerline rather than averaged across the outlet to avoid the spurious errors that are typically induced near mesh surfaces by virtue of vortex-ring singularities. By way of illustration, Fig. 12 compares the computed surface-vorticity centerline velocity to the one-dimensional flow estimation.

As expected, the potential solution based on the surface-vorticity formulation is found to be nearly one-dimensional. Nonetheless, aft of the finned section, approximately 4.44 m (175 in.) downstream, the surface-vorticity prediction leads to slightly higher Mach numbers than the one-dimensional prediction, thus confirming that the actual mass flow is more concentrated near the central axis.

To further verify the accuracy of the potential formulation, a streamtube analysis is performed. One may recall that the surface-vorticity solution is fundamentally isentropic, as there are no irreversibilities present in the irrotational and inviscid fields. Moreover, because the motion is tangent to every point along a streamline, no flow crosses a streamtube boundary. As such, the cross-sectional area and Mach number profile of a given streamtube must remain consistent with isentropic flow theory. A streamline near the centerline may be chosen along with discrete points along its length. Assuming axisymmetric flow, the streamline may be revolved around the central axis to create a streamtube with a finite circular cross section. Data near the central axis can then be used to avoid the known error near mesh surfaces and the noncircular streamtube cross sections in the finocyl section. At each of the discrete points defined along the streamline, the Mach number and streamtube's cross-sectional area can be recorded. These values can then be compared to Stodola's area–Mach relation from isentropic flow theory.⁹⁰ The latter is given by:

In the above, a reference point near the headwall is used, denoted by the subscript h, as there is no choke point in the computational model. After some algebra, a comparison between the streamtube analysis and the isentropic flow approximation can be carried out, thus leading to excellent agreement, as illustrated in Fig. 13. Unsurprisingly, the surface-vorticity flow remains isentropic throughout the chamber. Note that this particular verification step confirms that the mathematical solution to the equations of motion is consistent with the model's underlying assumptions; it is not intended to serve as a validation of the physical flow solution.

At this stage, one may proceed to compare the surface-vorticity's corrected solution to two computational models as well as the experimental pressure data acquired from RSRM tests. The two computational models are actually developed by Thiokol: the first one, SHARP, consists of a fully coupled compressible Navier–Stokes solver and the second, known as SCB02, refers to Thiokol's internal one-dimensional code.⁸⁹ The experimental pressure measurements are gathered from two RSRM qualification motors, namely, QM-7 and QM-8. In this comparison, the measurements from the qualification motors are rescaled to account for the difference in the propellants' initial temperature and burn rate. Forthwith, Mach number and pressure distributions from these computational and experimental data are showcased in Fig. 14 along the RSRM's central axis.

It may be instructive to note that, in these simulations, the rotational and compressibility correction are applied at the outlet of the finocyl, i.e., at approximately 4.44 m (175 in.) from the headwall. As a result, the velocity may be seen to markedly increase at the outlet of the finocyl. Particularly, the velocity increases by approximately 3% in view of the customized inlet condition. Moreover, one finds that the effect of compressibility correction is to elevate the velocity further by approximately 5%. The resulting comparisons, provided in Fig. 14, confirm that the irrotational solution agrees well with the other computational models in the finocyl section. This behavior could be partly attributed to the lack of rotational energy in the central bore of the finocyl. The surface-vorticity calculations also agree well with the Mach number profiles of the SHARP solver as well as the pressure profiles of SHARP, SCB02, and the QM-7 motor in the first half of the RSRM.

However, the surface-vorticity solution appears to lose accuracy in the aft half of the RSRM, particularly, at the locations of the radial slots separating the motor segments; the effects of these are observed in the SHARP and SCB02 data, specifically as the two sharp gradients located at approximately 17 and 25 m, respectively. This behavior may be further attributed to the slots having a much larger influence in a rotational than an irrotational flow. In a rotational field, the injected stream from the radial slots induces a rise in pressure just forward of the slot and a suction pressure just aft of the slot. A recirculation zone generally forms aft of the slot as a resultant of the suction pressure.⁹¹ The recirculating flow restricts the bulk flow, thus leading to an increase in its mass flow rate. In irrotational flow, the radial slots represent virtually thin strips along the burning surface with slightly increased mass addition. One can infer from Fig. 14 that the surface-vorticity solution along the centerline is affected only slightly by the presence of radial slots. Given that the rotational correction consists of a linear augmentation to the irrotational flow, the presence of slots is not adequately captured in the surface-vorticity results. Nonetheless, the slot between the first two motor segments happens to be much smaller than the other two slots; as such, it remains inconsequential to the centerline motion in both the irrotational and rotational solutions.

A possible avenue to account for the additional rotationality induced by the intersegmental gaps and inhibitors may be based on the concept of “vortex walls”; these virtual boundaries act as impermeable walls to the potential solution and can be prescribed by flow recirculation and vortices generated by the grain inhibitors. Such vortex walls have been effectively used in conjunction with a potential solver by Smith-Kent et al.⁸ in their modeling effort of the cavity formed around a submerged nozzle. The difficulty in utilizing this approach, however, may be associated with the degree of precision entailed in measuring the recirculation region along with the mass transfer model needed to represent the diffusive mixing between the recirculation zone and the bulk flow. For demonstrative purposes, an estimated recirculation zone boundary (RZB) is applied to the RSRM model. Its incorporation leads to markedly improved results, as shown in Fig. 14, using dash double-dotted lines. In this analysis, the forward half of the model remains an injecting surface boundary, that is, unaltered from the initial propellant geometry. The aft section, however, can be used to incorporate an RZB region, as shown in Fig. 15, namely, beginning at the middle motor slot, as a virtual boundary between recirculating fluid and the main flow. Presently, there is no model in the surface-vorticity solver to accurately predict the mass transfer rate across the RZB. For simplicity, the surface normal injection speed at the boundary is used to represent this estimate. Assuming that the recirculation zone is quasi-steady, the mass flux into the RZB domain from the burning surface can be set equal to the flow exiting the RZB and merging with the bulk flow as a result of diffusion. This arrangement enables the model to satisfy mass conservation across the RZB. Moreover, considering that the injection Mach number is of order $0.01−0.001$⁠, the error introduced by the assumption of a uniform mass transfer rate across the RZB may be neglected. The geometry of the RZB can be approximated rather straightforwardly as an offset boundary from the burning surface, as shown in Fig. 15. Note that the first portion of the RZB starts to develop from the second radial slot (⁠$z¯=16.8 m$⁠) and then spans the length of the third motor segment. The RZB geometry in the third motor segment runs parallel to the burning surface with a displacement thickness of $6.35$ cm. To terminate the RZB domain in the third motor segment, a spline curve may be applied between the offset boundary and the burning surface at the radial slot. It is apparent from the experimental data as well as the SHARP simulations that the rotationality introduced into the flow by the third radial slot (at $z¯=24.8 m$⁠) further restricts the bulk motion. Bearing this in mind, the RZB is regenerated for the fourth motor segment using a similar strategy and a larger displacement thickness of 11.9 cm. Since the aft section of the fourth motor segment exposes a relatively large taper, the adverse pressure gradient induced by this taper increases the thickness of the recirculation zone; this effect causes the bulk flow to experience a smaller taper angle than that prescribed by the actual propellant surface. Consequently, the RZB taper in the aft section can be set to a lower angle of $2.5°$ Note that the rotational effect of a submerged nozzle, had one been considered, would have been to further restrict the bulk flow at the outlet of the chamber. As seen in Fig. 14, the inclusion of RZB leads to closer agreement with the experimental measurements in the aft section of the RSRM. In future work, it is hoped that the determination of the RZB geometry will be guided by rigorous theoretical considerations even in the absence of experimental measurements. The same may be said of the mass transfer model across the RZB, whose formulation and refinement remain an open research question.

In returning to Fig. 14, it can be noted that, in the absence of an RZB correction, a sharp gradient may be detected in the surface-vorticity computations at the outlet of the motor; this occurs in the tapered section of the last motor segment. One may recall that the rotational correction, that is, incorporated here is based on analytical models that correspond to a small taper angle. The taper angle in the motor segment becomes approximately 10°, thus leading to more noticeable deviations. For such large taper angles, the accuracy of the present model deteriorates. Another contributing factor may be ascribed to the submerged nozzle and the resultant cavity between the grain and the nozzle not being modeled in this analysis. The cavity likely behaves similarly to the radial slots as fluid is injected from the cavity into the main port, where the injected stream restricts the upstream mass flow and drives the flow rate to slightly higher values. Again, such effects could be more easily characterized using rotational solvers and more intensive correction methods. However, as with any fast capturing conceptual and preliminary design tool, one must be willing to trade some degree of accuracy to gain a substantially higher speed of execution. For example, the surface-vorticity computations carried out in this study using a standard personal computer have required several minutes to half an hour depending on the model's complexity; in contrast, the use of full Navier–Stokes simulations on high speed clusters can be orders of magnitude more computationally intensive and, as such, of more limited use in supporting iterative design permutations.⁹² 

The objective of this work has been to demonstrate the viability of using surface panel solver approaches for internal flow analysis, with particular interest in solid rocket motor chamber configurations. By relying on the surface-vorticity approach with appropriate modifications, high quality instream solutions have been realized for limiting process configurations that possess analytical solutions. These include cylindrical and planar geometries with slight taper and headwall injection. In this study, one verifies that the solutions realized using the enhanced surface-vorticity solver agree well with the analytical models for all of the cases tested. The error near the burning surface is also anticipated due to the singularities present on the mesh surface, an inevitable characteristic of panel solvers. In fact, work is presently underway to develop a vortex-splitting scheme that aims at reducing this error. In the planar models, one observes a slight error around the midsection plane that seems to grow in magnitude as the length of the chamber is further increased. This is possibly due to the finite width of the models. To avoid edge effects in the spanwise direction, the width of the porous channel has been increased to at least six times its half height, as suggested by Terrill,⁸⁷ to the extent of reducing these errors substantially.

In the RSRM validation study, the potential solution has been first verified through both continuity-preserving and streamtube analyses. The comparisons have been rather favorable: the irrotational motion follows quite closely the one-dimensional approximations in a fundamentally isentropic manner. A more physical comparison has been carried out using experimental data and other computational models developed by Thiokol. Despite the complex internal environment of the RSRM, the surface-vorticity solver is shown to reproduce accurate Mach and pressure profiles with only two simple flow corrections. It is demonstrated through this effort that the surface-vorticity solution in the forward half of the RSRM matches exceptionally well the more rigorous SHARP solver and the pressure data of the QM-7 qualification motor. Nonetheless, the present approach does not predict the flow as robustly in the aft half, due in large part to recirculation caused by the radial slots. Similarly, the effects of the submerged nozzle are not incorporated into the present solution. Since the rotational correction uses a small taper angle approximation, it may not be as suitable for modeling motors with large aft motor tapers.

The next logical step in the internal ballistics analysis using this mid-fidelity approach is to couple grain regression with the enhanced surface-vorticity technique to produce thrust-time profiles. The regressive evolution of the burning surface could also be useful in erosive burning analyses to investigate its effect on the flowfield. Other internal flow patterns, such as those arising in hybrid rocket motors, could also be explored. The irrotational approach might prove useful in hybrid rocket applications where a limited-diffusion flame is present.

In a different context, acoustic studies, or other perturbation-based methods, could use the steady-state solution obtained from a surface-vorticity approach as a leading-order approximation or a basis upon which transient behavior could be superimposed. Many analytical studies based on potential flow theory could be extended to more complex geometries using the approach presented here. The development of internal ballistics and enhanced burn rate modeling capabilities that are implementable at the conceptual and preliminary stages of design strongly support the ongoing multidisciplinary design analysis and optimization (MDO) pipelines.

In hindsight, the irrotational flow assumption stands as the most significant limitation of the present solver in the treatment of solid rocket internal flowfields. As seen through the RSRM analysis, rotational effects can have a considerable impact on the internal flowfield. In this study, rotationality is incorporated into the solution along the centerline using a simple correction that does not encompass the entire chamber; it also omits the effect of the radial slots. Other methods may be explored to approximate, or fully integrate, rotational flow effects into potential flow solutions that are based on surface meshes. A possible correction presently being investigated leverages the concept of a recirculation zone boundary (RZB). In this work, the surface-vorticity predictions that incorporate an estimated RZB correction are shown to clearly outperform their strictly irrotational counterparts for every RSRM test case. To further improve the procedure needed to integrate an RZB region into a panel code, it will be helpful to accurately capture the shape of the recirculation zone. The same may be said of a diffusive mass transfer model from the RZB to the bulk fluid motion whose determination remains an open subject of investigation.

The authors report no conflict of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.