Open AccessRegular Article

Viscous Mean Flow Approximations for Porous Tubes with Radially Regressing Walls

Tony Saad

University of Utah, Salt Lake City, Utah 84112

*Assistant Professor, Department of Chemical Engineering.

Search for more papers by this author

and

Joseph Majdalani

Auburn University, Auburn, AL 36849

†Professor and Francis Chair of Excellence, Department of Aerospace Engineering. Associate Fellow AIAA.

Search for more papers by this author

Published Online:17 Jul 2017https://doi.org/10.2514/1.J055949

Abstract

The mean gaseous motion in solid rocket motors has been traditionally described using an inviscid solution in a porous tube of fixed radius and uniform wall injection. This model, usually referred to as the Taylor–Culick profile, consists of a rotational solution that captures the bulk gaseous motion in a frictionless rocket chamber. In practice, however, the port radius increases as the propellant burns, thus leading to time-dependent effects on the mean flow. This work considers the related problem in the context of viscous motion in a porous tube and allows the radius to be time dependent. By implementing a similarity transformation in space and time, the incompressible Navier–Stokes equations are first reduced to a nonlinear fourth-order ordinary differential equation with four boundary conditions. This equation is then solved both numerically and asymptotically, using the injection Reynolds number $R e$ and the dimensionless wall regression ratio $α$ as primary and secondary perturbation parameters. In this manner, closed-form analytical solutions are obtained for both large and small Reynolds number with small-to-moderate $α$ . The resulting approximations are then compared with the numerical solution obtained for an equivalent third-order ordinary differential equation in which both shooting and the irregular limit that affects the fourth-order formulation are circumvented. This code is found to be capable of producing the stable solutions for this problem over a wide range of Reynolds numbers and wall regression ratios.

Nomenclature
$a$	time-dependent radius
$\dot{a}$	wall regression speed
$b$	scaling factor
$F$	characteristic mean flow function
$G$	transformed mean flow function
$K$	characteristic integration constant, $K (R e)$
$n_{\max}$	total number of points in the integration domain
$n_{p}$	point where $G^{'} (n_{p}) = 0$
$n_{Taylor}$	number of points for which the Taylor series expansion is applied
$R e$	Reynolds number
$r$	radial coordinate
$S$	surface
$t$	time
$U_{m}$	mean axial velocity
$u$	velocity vector
$V$	absolute injection velocity as seen by a stationary observer
$V$	volume
$V_{w}$	injection velocity with respect to the moving wall
$z$	axial coordinate
$α$	wall expansion ratio
$β$	scaled expansion ratio, $α / b$
$γ_{i}$	$i$ th coefficient in the Taylor series expansion
$ε$	small real number used in numerical integration
$ϵ$	small perturbation parameter
$η$	radial transformation variable, $(1 / 2) r^{2}$
$λ$	scaling factor, $2 / R e$
$ν$	kinematic viscosity
$ξ$	transformation variable, $b η$
$ξ_{\max}$	domain integration size
$ξ_{p}$	coordinate location corresponding to $n_{p}$
$ρ$	density
$ψ$	Stokes stream function
$Ω$	vorticity vector

I. Introduction

Viscous motion in cylindrical chambers with sidewall injection is of interest in a variety of applications, including mean flow modeling of solid [1–3] and hybrid rockets [4,5], sweat cooling [6,7], boundary-layer control [8–10], peristaltic pumping [11,12], gaseous diffusion, and isotope separation [13–15]. It is the latter group of studies by Berman [13,14] that has actually provided the impetus to develop the first similarity transformation of the Navier–Stokes equations into a fourth-order nonlinear ordinary differential equation (ODE). The resulting ODE was later solved asymptotically over different ranges of the wall injection Reynolds number. Most notable examples include those by Yuan and Finkelstein [16], Terrill and Thomas [17], Terrill [18,19], Skalak and Wang [20], and others.

In the propulsion area, an inviscid rotational counterpart known as the Taylor [21], Culick [22], or Taylor–Culick solution occupies the central stage in modeling the bulk gaseous motion inside solid rocket motors (SRMs). This may be owed to its association with several studies involving hydrodynamic instability [23 –28], acoustic instability [29 –35], wave propagation [36 –39], particle-mean flow interactions [40], and rocket performance measurements [41–43]. The Taylor–Culick solution was originally verified to be an adequate representation of the expected flowfield in SRMs both numerically by Sabnis et al. [44] and experimentally by Dunlap et al. [45,46], thereby confirming its character in a nonreactive chamber environment. It was extended by Majdalani and Akiki [4] to include effects of viscosity and headwall injection, by Saad et al. [42] and Sams et al. [43] to account for wall taper, by Kurdyumov [47] to capture effects of irregular cross sections, and by Majdalani and Saad [48] to allow for arbitrary headwall injection. Then, using variational calculus and the Lagrangian optimization principle, Saad and Majdalani [49] uncovered a continuous spectrum of Taylor-like solutions exhibiting increasing kinetic energy signatures, while ranging from the traditional Culick profile down to its predecessor, the irrotational mean flow known as the Hart–McClure profile [50,51]. This potential mean flow preceded the use of Culick’s model in several fundamental investigations of combustion instability [52–54]. However, these studies have not been concerned with the time dependence that is inherent to the grain surface of a burning propellant.

In a practical analysis, the radius of a rocket motor is allowed to increase during wall regression, and the corresponding problem gives rise to a partial differential equation (PDE) in lieu of an ODE. Although the problem in a tube with expanding or contracting surfaces was treated by Uchida and Aoki [12] in the late 1970s, the effects of an injecting sidewall were incorporated more than a decade later by Goto and Uchida [55], in the context of a pulsating porous tube and by Majdalani et al. [1,2] in the context of a rocket chamber. In their work, the effect of wall regression was prescribed by a dimensionless wall expansion ratio, $α$ , written as a viscous Reynolds number based on the radial regression speed of the sidewall. In the interim, Dauenhauer and Majdalani [56–58] formulated the corresponding equation and numerical solution for the planar flow analog. The analytical solutions promoted through these efforts would later receive attention in follow-up studies aimed at devising asymptotic approximations over various ranges of the control parameters or at reconstructing the solution using alternative techniques such as the Lie-group theory [59,60] or the Homotopy-Analysis Method [61–63].

In this article, we revisit the problem of viscous motion in a porous tube and allow the radius to transiently expand or contract [1,2]. As before, we employ a dual similarity transformation in space and time to reduce the Navier–Stokes equations into a nonlinear fourth-order ODE that can be solved both asymptotically and numerically. After reconstructing the analytical solution of Majdalani et al. [1,2] for large injection, we employ two perturbation parameters, the injection Reynolds number $R e$ and the wall expansion ratio $α$ to formulate an asymptotic series for the case of small Reynolds number and small-to-moderate $α$ . We then introduce an efficient numerical approach that overcomes the singularity encountered in this model, specifically in the form of an intrinsically satisfied limit at the origin. To make headway, an equivalent third-order ODE is presented and solved using a technique that may be attributed to Terrill and Thomas [17]. Subsequent numerical predictions are obtained directly, with no shooting or iteration, and then compared with our analytical approximations. In this manner, the robustness of the numerical algorithm is used to demonstrate that the present approach can capture the stable solutions for this problem over a wide range of Reynolds number and $α$ . Although the case of large Reynolds number is connected with modeling solid rocket motors, the small injection case is relevant to other problems involving sweat cooling, filtration, peristaltic pumping, and boundary-layer control.

The paper is organized as follows. First, we consider the mathematical model developed by Majdalani et al. [1,2] and reaffirm its validity. Second, we discuss the numerical technique that appears to be most practical. Third, we produce regular perturbation approximations for the large and small injection cases, the latter of which being are presented for the first time. Finally, we compare asymptotics and numerics over a common range of the control parameters.

II. Mathematical Model

The cylindrical propellant grain of a solid rocket motor is modeled as a long tube with one end closed at the headwall and the other end open. The circumference of the wall is assumed to be permeable so as to simulate propellant burning and normal gas injection. Furthermore, the wet area of the tube is allowed to radially expand at a speed equal to $\dot{a}$ . For this to occur, the mathematical model requires that the headwall area stretches accordingly to accommodate the expansion of the tube. This behavior is shown in Fig. 1, where an axisymmetric coordinate system is defined. By assuming an incompressible mean flow, the vorticity transport equation is written as

\frac{\partial Ω^{*}}{\partial t} - \nabla^{*} \times (u^{*} \times Ω^{*}) = ν \nabla^{*} \times \nabla^{* 2} u^{*}

(1)

Fig. 1
Schematic of a cylindrical chamber used to illustrate sidewall injection and wall regression, as well as the control volume used to evaluate the average axial velocity.

A. Boundary Conditions

These can be organized as follows:

r^{*} = a (t) : {\begin{matrix} u_{r}^{*} = - V & sidewall injection \\ u_{z}^{*} = 0 & no slip \end{matrix}

(2)

r^{*} = 0 : {\begin{matrix} u_{r}^{*} = 0 & no flow across the centerline \\ \frac{\partial u_{z}^{*}}{\partial r^{*}} = 0 & axial velocity symmetry \end{matrix}

(3)

z^{*} = 0 : u_{z}^{*} = 0 inert headwall

(4)

Here,

V

is the absolute fluid injection velocity at the wall. It is defined as the fluid velocity seen by an observer in a reference frame located outside of the tube. If

V_{w}

is used to denote the fluid velocity with respect to the wall, then the following relation holds:

V = V_{w} - \dot{a}

(5)

Evidently, when the walls are stationary, the absolute and relative velocities become equal. But when the wall undergoes inward contraction, then $\dot{a} < 0$ and the fluid will appear as though it is being injected at a larger speed than $V_{w}$ . Conversely, if the wall regresses outwardly, then $\dot{a} > 0$ and the fluid injection velocity will appear to be smaller than $V_{w}$ .

Note that, at $z^{*} = 0$ , we do not enforce any condition on the radial velocity profile. This slip condition on the radial velocity is what actually allows us to solve this problem by assuming a linear relation between the stream function and the axial coordinate $z^{*}$ [see Eq. (14)]. In reality, a small, rather negligible viscous layer develops at the headwall, whose analysis is described by Chedevergne et al. [24].

B. Similarity in Space

By inspection, one expects the axial velocity to vary linearly in $z^{*}$ . To show this, we follow Majdalani et al. [1–3] by considering the control volume delineated by the dashed lines in Fig. 1. Through a mass balance on the inlet and outlet surfaces, we require

\frac{\partial V}{\partial t} + \iint_{S_{exit}} u^{*} \cdot n d S - 2 π V_{w} a z^{*} = 0

(6)

The term containing the integral in Eq. (6) may be connected to the average axial velocity $U_{m}$ inside the volume using

U_{m} = \frac{1}{π a^{2}} \iint_{S_{exit}} u^{*} \cdot n d S

(7)

Its substitution into Eq. (6) yields

\frac{\partial V}{\partial t} + π a^{2} U_{m} - 2 π V_{w} a z^{*} = 0

(8)

and so, by using

V = π a^{2} z^{*}

, Eq. (8) reduces to

2 π a \dot{a} z^{*} + π a^{2} U_{m} - 2 π V_{w} a z^{*} = 0

(9)

Thus, based on mass conservation alone, we must have

U_{m} = 2 (V_{w} - \dot{a}) \frac{z^{*}}{a} = 2 V \frac{z^{*}}{a}

(10)

which must be equated to the area integral for the average velocity

U_{m} (z^{*}, t) = \frac{1}{π a^{2}} \int_{0}^{a} u_{z}^{*} (r^{*}, z^{*}, t) r^{*} d r^{*} = 2 V \frac{z^{*}}{a}

(11)

which proves that the only admissible form of the axial velocity is a linear relation of the type

u_{z}^{*} (r^{*}, z^{*}, t) = z^{*} f (r^{*}, t)

(12)

C. Vorticity Transport Equation

We introduce the Stokes stream function using

u_{r}^{*} = - \frac{1}{r^{*}} \frac{\partial ψ^{*}}{\partial z^{*}}; u_{z}^{*} = \frac{1}{r^{*}} \frac{\partial ψ^{*}}{\partial r^{*}}

(13)

Then, based on Eq. (12), the stream function may be written as

ψ^{*} = ν z^{*} F (r, t); r = \frac{r^{*}}{a}

(14)

where

ν

denotes the kinematic viscosity for a dimensionless

F (r, t)

. The radial and axial velocities can now be expressed in terms of

F

. This operation yields

u_{r}^{*} = - \frac{ν}{a r} F, u_{z}^{*} = \frac{ν z^{*}}{a^{2} r} \frac{\partial F}{\partial r}

(15)

Because $u_{r}^{*}$ remains independent of the axial coordinate, the vorticity displays a single nonzero component in the tangential direction, viz.,

Ω^{*} = Ω_{θ}^{*} = \frac{\partial u_{r}^{*}}{\partial z^{*}} - \frac{\partial u_{z}^{*}}{\partial r^{*}} = - \frac{\partial u_{z}^{*}}{\partial r^{*}}

(16)

Upon substitution of Eq. (16) into Eq. (1), one recovers

\frac{\partial}{\partial r^{*}} (\frac{\partial u_{z}^{*}}{\partial t}) + \frac{\partial}{\partial r^{*}} (u_{r}^{*} \frac{\partial u_{z}^{*}}{\partial r^{*}}) + \frac{\partial}{\partial r^{*}} (u_{z}^{*} \frac{\partial u_{z}^{*}}{\partial z^{*}}) - ν \frac{\partial}{\partial r^{*}} [\frac{1}{r^{*}} \frac{\partial}{\partial r^{*}} (r^{*} \frac{\partial u_{z}^{*}}{\partial r^{*}})] = 0

(17)

By paying careful attention to partial differentiation, Eq. (17) is reduced to

{- \frac{a^{2}}{ν} \frac{F_{r t}}{r} + \frac{a \dot{a}}{ν} r {(\frac{F_{r}}{r})}_{r} + 2 \frac{a \dot{a}}{ν} \frac{F_{r}}{r} + \frac{F}{r} {(\frac{F_{r}}{r})}_{r} - {(\frac{F_{r}}{r})}^{2} + \frac{1}{r} {(\frac{F_{r}}{r})}_{r} {(\frac{F_{r}}{r})}_{r r}}_{r} = 0

(18)

where the subscripts

r

and

t

denote partial differentiation with respect to

r

and

t

, respectively. At this juncture, we define the wall expansion ratio

α (t) \equiv \frac{a \dot{a}}{ν}

(19)

and convert Eq. (18) into

{{(\frac{F_{r}}{r})}_{r r} + (\frac{1}{r} + \frac{F}{r} + α r) {(\frac{F_{r}}{r})}_{r} - (\frac{F_{r}}{r} - 2 α) \frac{F_{r}}{r} - \frac{a^{2}}{ν} \frac{F_{r t}}{r}}_{r} = 0

(20)

The outcome can be readily integrated with respect to $r$ and rearranged into

{(\frac{F_{r}}{r})}_{r r} + (\frac{1}{r} + \frac{F}{r} + α r) {(\frac{F_{r}}{r})}_{r} - (\frac{F_{r}}{r} - 2 α) \frac{F_{r}}{r} - \frac{a^{2}}{ν} \frac{F_{r t}}{r} = K

(21)

The formidable PDE that we arrive at embodies the physics of viscous motion in a porous pipe with retracting walls. It is related to work by Uchida and Aoki [12], who studied a similar problem in the context of pipe flow with expanding or contracting impermeable walls. Goto and Uchida [55], Dauenhauer and Majdalani [56–58], and Majdalani et al. [1,2] later extended Uchida and Aoki [12] to permit fluid injection or suction along the moving sidewall. Their analysis was carried out for both planar [56–58] and cylindrical configurations [1,2]. Majdalani et al. [1–3] also derived asymptotic solutions for this equation over a practical range of Reynolds numbers and expansion ratios.

D. Similarity in Time

At first glance, Eq. (21) may seem intractable using any of the standard analytical techniques. This excludes the Homotopy Analysis Method, which has been shown to handle highly nonlinear equations quite admirably [61,64–66]. By using a practical hypothesis that may be traced back to Uchida and Aoki [12], Eq. (21) may be reduced to a third-order nonlinear ODE. As shown by Majdalani et al. [1], this is accomplished by setting

α (t) = constant; F (r, t) \mapsto F [r, α (t)]

(22)

The constancy of $α$ ensures that

F_{r t} = \frac{\partial F_{r}}{\partial α} \frac{d α}{d t} = 0

(23)

One can then deduce the required wall regression speed for which $α$ becomes a constant. According to this model

α = \frac{a \dot{a}}{ν} = \frac{a_{0} {\dot{a}}_{0}}{ν}

(24)

where

a_{0}

and

{\dot{a}}_{0}

correspond to the initial radius and regression speed, respectively. Upon integration, we get

a (t) = a_{0} \sqrt{1 + 2 \frac{ν α}{a_{0}^{2}} t}

(25)

Note that this model does not suppress time dependence, but rather embeds it within the solution implicitly. This is realized by specifying the regression speed in such a way that $α (t)$ remains constant.

Before substituting the time similarity conditions into Eq. (21), we find it useful to introduce the following normalizations:

F \mapsto \frac{F}{R e}; η \equiv \frac{1}{2} r^{2}

(26)

where

R e

is the sidewall injection Reynolds number based on the absolute velocity of the fluid:

R e \equiv \frac{V a}{ν} = \frac{V_{w} a}{ν} - α

(27)

In this setting, the Reynolds number remains a function of time. To make further headway, we assume that $R e = constant$ or

\frac{V_{w} a}{ν} = \frac{V_{w 0} a_{0}}{ν}

(28)

This requirement enables us to extract the proper form of $V_{w}$ in our model, namely,

V_{w} = \frac{V_{w 0} a_{0}}{a} = V_{w 0} {(1 + 2 \frac{ν α}{a_{0}^{2}} t)}^{- (1 / 2)}

(29)

Finally, backward substitution of Eq. (26) into Eq. (21) yields

η F^{'''} + F^{''} + \frac{1}{2} R e (F F^{''} - F^{' 2}) + α (η F^{''} + F^{'}) = K (R e)

(30)

where primes denote differentiation with respect to

η

. To eliminate the constant

K (R e)

, one differentiates Eq. (30) with respect to

η

and writes

η F^{''''} + α (η F^{'''} + 2 F^{''}) + \frac{1}{2} R e (F F^{'''} - F^{'} F^{''}) + 2 F^{'''} = 0

(31)

Finally, the boundary conditions may be assimilated into the following set

\frac{d F (1 / 2)}{d η} = 0

(32a)

F (\frac{1}{2}) = 1

(32b)

F (0) = 0

(32c)

\lim_{η \to 0} \sqrt{2 η} \frac{d^{2} F}{d η^{2}} = 0

(32d)

Note that the constant $K (R e)$ may be specified by evaluating Eq. (30) at any location inside the domain, such as the centerline. At $r = 0$ , one gets

K (R e) = F^{''} (0) + \frac{1}{2} R e F^{' 2} (0) + α F^{'} (0)

(33)

Either of the differential equations given by Eqs. (30) and (31) may be used to solve for the mean flow function $F (η)$ . The choice depends on the type of solution sought. For example, we elect to use Eq. (30) for the numerical solution because it requires less memory storage and fewer Runge–Kutta integrations. On the other hand, Eq. (31) will be used to initiate the regular perturbation expansion because of the implicit nonlinearity that the constant $K (R e)$ embodies.

III. Numerical Technique

Being nonlinear ODEs, Eqs. (30) and (31) may be solved directly through a Runge–Kutta integration routine. However, a close inspection of the boundary conditions as well as the governing equations reveals several difficulties. First, the explicit presence of $η$ in these equations requires a careful treatment by virtue of the singularity at the origin. This matter may be settled through a Taylor series expansion near the centerline. Second, the boundary conditions given by Eq. (32) result in a boundary value problem that requires a double infinity of integrations [17]. This problem may be overcome via a shooting method based on a Newton–Raphson root finding algorithm, such as the one proposed by Dauenhauer and Majdalani [58]; however, the fourth boundary condition in Eq. (32d) is difficult to implement. The reason is that, as long as $F^{''} (0)$ remains finite, Eq. (32d) will be automatically satisfied. Indeed, one expects $F^{''} (0)$ to be finite because it corresponds to the vorticity at the centerline. In short, we find Eq. (32d) to be impractical for numerical treatment.

The aforementioned issues may be overcome if one uses a transformation that is nearly identical to that employed by Terrill and Thomas [17]. These researchers encountered difficulties in an analogous problem involving flow through a porous tube. In our model, however, one needs to incorporate the effect of $α$ . By allowing both the Reynolds number and $α$ to be determined a posteriori (i.e., at the end of integration), a shooting method will no longer be required and the numerical solution may be arrived at directly, using a single Runge–Kutta integration. To illustrate this procedure, we begin by introducing the dual transformation

F = λ G (ξ); ξ = b η

(34)

where both

λ

and

b

are scaling factors that may be determined once the numerical integration is complete. Upon substitution into Eq. (30), we recover a third-order ODE in

G (ξ)

, namely,

ξ G^{'''} + G^{''} + \frac{1}{2} λ R e (G G^{''} - G^{' 2}) + \frac{α}{b} (ξ G^{''} + G^{'}) = \frac{K}{λ b^{2}} = K_{1}

(35)

where the primes denote differentiation with respect to

ξ

. If we further set

λ = \frac{2}{R e} and β = \frac{α}{b}

(36)

then Eq. (35) reduces to the convenient form

ξ G^{'''} + G^{''} + G G^{''} - G^{' 2} + β (ξ G^{''} + G^{'}) = K_{1}

(37)

Finally, by applying this procedure to the conditions encapsulated in Eq. (32), they transform into

\frac{d G (b / 2)}{d ξ} = 0

(38a)

G (\frac{1}{2} b) = \frac{1}{λ} = \frac{R e}{2}

(38b)

G (0) = 0

(38c)

\lim_{ξ \to 0} \sqrt{ξ} \frac{d^{2} G}{d ξ^{2}} = 0

(38d)

By initializing the solution with arbitrary values for $G^{'} (0)$ , $G^{''} (0)$ , and $β$ , the integration may be carried out until $G^{'} (b / 2) = 0$ . In practice, to produce a desired Reynolds number and wall expansion ratio, these are chosen based on the asymptotic solutions (to be derived in the next section). The point at which $G^{'}$ vanishes enables us to explicitly deduce $b$ with no need for iteration. Subsequently, the Reynolds number and $α$ may be calculated from Eq. (38b) as $R e = 2 G (b / 2)$ and $α = β b$ . The constant $K_{1}$ can also be extracted from the initial guesses by putting

K_{1} = G^{''} (0) - G^{'} {(0)}^{2} + β G^{'} (0)

(39)

Finally, near centerline ( $ξ \to 0$ ), a Taylor series expansion of the governing equation is required. Using

G (ξ) = \sum_{i = 0} γ_{i} ξ^{i}

(40)

substitution into Eq. (37) yields

γ_{0} = G (0) = 0; γ_{1} = G^{'} (0); γ_{2} = \frac{1}{2} G^{''} (0)

(41)

Note that these coefficients are known because the values of the first and second derivatives at the centerline are used to seed the numerical solution. Then, by substituting Eq. (40) into Eq. (30), coefficients of the same order may be equated to the extent of establishing a recurrence formula for the coefficients of the Taylor series expansion; we get

(n + 2) {(n + 1)}^{2} γ_{n + 2} = \sum_{j = 0}^{n} (j + 1) (n - j + 1) γ_{j + 1} γ_{n - j + 1} - \sum_{j = 0}^{n} (j + 1) (j + 2) γ_{j + 2} γ_{n - j} - β {(n + 1)}^{2} γ_{n + 1}; n \geq 1

(42)

Algorithm 1 lists the necessary steps for our numerical implementation. Results of the numerical simulation will be presented in subsequent sections, where they will be compared with the asymptotic approximations.

IV. Asymptotic Treatment

A. Large Injection

In this case, which is appropriate for modeling gaseous injection in solid rocket motors, the Reynolds number is sufficiently large to warrant a regular perturbation in the inverse of the Reynolds number. We follow Majdalani et al. [1] and write

F = F_{0} + ϵ F_{1} + O (ϵ^{2}); ϵ \equiv \frac{1}{R e}

(43)

Upon substitution into Eq. (31), we recover the ODEs for the leading and first-order solutions. These are

O (0) : F_{0} F_{0}^{'''} - F_{0}^{'} F_{0}^{''} = 0

(44a)

O (1) : F_{0} F_{1}^{'''} + F_{1} F_{0}^{'''} - F_{0}^{'} F_{1}^{''} - F_{1}^{'} F_{0}^{''} + 2 η F_{0}^{''''} + (2 α η + 4) F_{0}^{'''} + 4 α F_{0}^{''} = 0

(44b)

Note that we limit our analysis to the leading and first-order approximations partly due to the excellent agreement that we later show between the present solution and numerical calculations.

1. Leading-Order Solution

At leading order, we have

F_{0} F_{0}^{'''} - F_{0}^{'} F_{0}^{''} = 0

(45)

This can be solved by inspection to obtain

F_{0} = \sin (π η) = \sin θ; θ = π η

(46)

Equation (46) corresponds to the Taylor–Culick mean flow with no wall regression [22].

2. First-Order Solution

By substituting Eq. (46) into Eq. (44b), the ODE for the first-order term may be determined:

\sin θ F_{1}^{'''} - \cos θ F_{1} - \cos θ F_{1}^{''} + \sin θ F_{1}^{'} = - 2 ϑ \sin θ + (2 π^{- 1} α θ + 4) \cos θ + 4 π^{- 1} α \sin θ

(47)

The solution of Eq. (47) requires the identification of a homogeneous part $F_{1 h}$ . This can be obtained by solving

\sin θ F_{1 h}^{'''} - \cos θ F_{1 h} - \cos θ F_{1 h}^{''} + \sin θ F_{1 h}^{'} = 0

(48)

for which a general homogeneous solution can be defined as

F_{1 h} = C (θ) \cos θ

(49)

Upon substitution of Eq. (49) into Eq. (48), one obtains the differential equation governing $C (θ)$ , namely,

C^{'''} \sin θ \cos θ - 2 C^{''} \sin^{2} θ - C^{''} = 0

(50)

The solution of Eq. (50) may be managed through division by $C^{''} \sin θ \cos θ$ before integration. This operation yields

C (θ) = \frac{1}{2} K_{0} \tan θ + K_{1} θ + K_{2}

(51)

Finally, the general homogeneous solution of Eq. (48) emerges as

F_{1 h} = C (θ) \cos θ = K_{0} \sin θ + K_{1} θ \cos θ + K_{2} \cos θ

(52)

To obtain the total solution for Eq. (47), we use variation of parameters and rewrite

F_{1} = K_{0} (θ) \sin θ + K_{1} (θ) θ \cos θ + K_{2} (θ) \cos θ

(53)

then, upon differentiation, we get

F_{1}^{'} = K_{0} \cos θ + K_{1} (\cos θ - θ \sin θ) - K_{2} \sin θ

(54a)

F_{1}^{''} = - K_{0} \sin θ + K_{1} (- 2 \sin θ - θ \cos θ) - K_{2} \cos θ

(54b)

F_{1}^{'''} = - 2 K_{1^{'}} \sin θ - K_{0} \cos θ + K_{1} (- 2 \cos θ - \cos θ + θ \sin θ) + K_{2} \sin θ

(54c)

along with the auxiliary conditions

K_{0}^{'} \sin θ + K_{1}^{'} θ \cos θ + K_{2}^{'} \cos θ = 0

(55a)

K_{0}^{'} \cos θ + K_{1}^{'} (\cos θ - θ \sin θ) - K_{2}^{'} \sin θ = 0

(55b)

Inserting Eq. (54) into Eq. (47) leads to

- 2 K_{1}^{'} \sin θ = (2 π^{- 1} α θ + 4) \cot θ + 4 π^{- 1} α - 2 θ

(56)

At this juncture, one can solve Eq. (56) in conjunction with Eq. (55) simultaneously for $K_{0}^{'}$ , $K_{1}^{'}$ , and $K_{2^{'}}$ . The resulting linear system of equations produces the following solution

K_{0} = \frac{α}{π} (\cos θ + 3 \ln \tan \frac{1}{2} θ - θ \sin θ - θ \csc θ) - θ \cos θ - \sin θ - 2 \csc θ - S (θ) + C_{0}

(57)

K_{1} = \frac{α}{π} (θ \csc θ - 3 \ln \tan \frac{1}{2} θ) + 2 \csc θ + S (θ) + C_{1}

(58)

K_{2} = \frac{α}{π} [3 S (θ) - \sin θ - θ^{2} \csc θ - θ \cos θ] - \cos θ + θ \sin θ - 2 θ \csc θ - S_{1} (θ) + C_{2}

(59)

where

S (θ) = \int_{0}^{θ} ϕ \csc ϕ d ϕ; S_{1} (θ) = \int_{0}^{θ} ϕ^{2} \csc ϕ d ϕ

(60)

and

C_{0}

C_{1}

, and

C_{2}

are constants that can be determined by invoking the boundary conditions in Eq. (32). The first-order solution given by Eq. (53) is now at hand, specifically,

F_{1} = \frac{α}{π} [3 (\sin θ - θ \cos θ) \ln (\tan \frac{θ}{2}) - 2 θ] + \frac{3}{π} α \cos θ S (θ) - 3 - S_{1} (θ) \cos θ + (θ \cos θ - \sin θ) S (θ) + C_{0} \sin θ + C_{1} θ \cos θ + C_{2} \cos θ

(61)

What is left is the application of the boundary conditions, and these require that

F_{1} (\frac{1}{2} π) = 0 or C_{0} = 3 + α + S (\frac{1}{2} π)

(62)

in conjunction with

\lim_{θ \to 0} F_{1} (θ) = 0; or C_{2} = 3

(63)

and finally

\frac{d F_{1} (1 / 2)}{d η} = π \frac{d F_{1} (π / 2)}{d θ} = 0

(64)

from which we extract

C_{1} = - 1 - \frac{6}{π} + \frac{2}{π^{2}} α + \frac{2}{π} S_{1} (\frac{1}{2} π) - S (\frac{1}{2} π) (\frac{6}{π^{2}} α + 1)

(65)

3. Verification

The total solution may now be assembled into

F = F_{0} + ϵ F_{1} + O (ϵ^{2}) = \sin θ + ϵ {\frac{α}{π} [3 (\sin θ - θ \cos θ) \ln (\tan \frac{θ}{2}) - 2 θ] + \frac{3}{π} α \cos θ S (θ) - 3 - S_{1} (θ) \cos θ + (θ \cos θ - \sin θ) S (θ) + C_{0} \sin θ + C_{1} θ \cos θ + C_{2} \cos θ} + O (ϵ^{2})

(66)

This result is identical to the one obtained by Majdalani et al. [1,2], where corresponding flow attributes such as velocity, pressure, and vorticity are described. In this study, we limit ourselves to the detailed mathematical derivation and verification using a robust numerical solution. To this end, Eq. (66) is compared with results from the Runge–Kutta integration for Reynolds numbers of 100, 500, and 1000 taken at different values of the wall expansion ratio, namely, at $α = 10$ , 20, $- 10$ , and $- 20$ .

Starting with the mean flow function $F$ , corresponding graphs are shown in Figs. 2 and 3. The first noticeable feature is the accuracy of the solution for large Reynolds number regardless of the wall expansion ratio. In fact, even for the relatively small value of $R e = 100$ , the curves appear to be visually indistinguishable unless one uses magnification on certain areas of the graphs. The other noteworthy feature is the effect of decreasing the wall expansion ratio. Specifically, the curves become increasingly more difficult to separate as one scans the plots sequentially, starting with the largest value of $α$ in Figs. 2b and 3b.

Fig. 2
Comparison between analytical and numerical solutions for the mean flow function $F$ using a) $α \approx 10$ and b) $α \approx 20$ . Curves are shown for $R e \approx 100$ (black), 500 (red), and 1000 (blue). Everywhere, inset frames are used with symbols for added clarity.

Fig. 3
Comparison between analytical and numerical solutions for the mean flow function $F$ using a) $α \approx - 10$ and b) $α \approx - 20$ . Curves are shown for $R e \approx 100$ (black), 500 (red), and 1000 (blue).

The second set of comparisons is provided in Figs. 4 and 5 for $F^{'}$ . We recall that the derivative of the mean flow function is directly connected to the axial velocity and, as such, the importance of arriving at a reliable analytical representation for $F^{'} (η)$ cannot be underrated. In this case, similar trends are observed as those corresponding to $F$ . The most prominent feature in these graphs may be the effect of $α$ on the centerline velocity. As the wall expansion ratio is decreased, so is the centerline velocity. Altogether, the curves are practically inseparable, a satisfactory observation given that our solution is of first-order accuracy.

Fig. 4
Comparison between analytical and numerical solutions for $F^{'}$ using a) $α \approx 10$ and b) $α \approx 20$ . Curves are shown for $R e \approx 100$ (black), 500 (red), and 1000 (blue).

Fig. 5
Comparison between analytical and numerical solutions for $F^{'}$ using a) $α \approx - 10$ and b) $α \approx - 20$ . Curves are shown for $R e \approx 100$ (black), 500 (red), and 1000 (blue).

Additional confirmatory comparisons between numerics and asymptotics are provided in Tables 1 and 2 where $F$ , $F^{'}$ , and $F^{''}$ are listed side-by-side along with their numerical predictions at three injection Reynolds numbers of 100, 500, and 1000, and either $α = 10$ (Table 1) or $- 10$ (Table 2). Corresponding values are provided at five characteristic radial positions corresponding to $η = 0, 0.125, 0.25, 0.375$ , and 0.5. Irrespective of whether the injecting surface is expanding or contracting, it may be seen in both tables that the accuracy in $F, F^{'}$ , and $F^{''}$ for the $R e = 1000$ case extends to 4, 3, and 2 significant digits, respectively. A slight degradation in precision occurs as the Reynolds number is lowered first to $R e = 500$ and then to 100. The overall agreement, however, continues to fall approximately within the universally accepted engineering tolerance of 5%.

Given the error entailed in the model itself, the first-order approximation presented here appears to be sufficiently adequate, especially in the modeling of SRM flowfields for which the Reynolds number can be quite large [24]. There may be situations, however, where higher accuracy is desired at intermediate ranges of Reynolds number, for example, between 1 and 100. For such cases, higher-order approximations will be required and these could be obtained by resorting to other techniques, such as Lie-group theory [59,60], Adomian decomposition [64], or Homotopy Analysis Method (HAM) [61].

B. Small Injection and Suction with Low Regression

In this case, both the Reynolds number and the expansion ratio are low. This necessitates rewriting Eq. (31) as

2 η F^{''''} + α (2 η F^{'''} + 4 F^{''}) + R e (F F^{'''} - F^{'} F^{''}) + 4 F^{'''} = 0

(67)

To solve Eq. (67) using a perturbation expansion, one needs to expand in both small parameters, Reynolds number and $α$ . This is done by setting

F = F_{0} + R e F_{1} + O (R e^{2})

(68)

F_{0} = F_{00} + α F_{01} + O (α^{2})

(69)

F_{1} = F_{10} + α F_{11} + O (α^{2})

(70)