Asymptotically based self-similarity solution of the Navier–Stokes equations for a porous tube with a non-circular cross-section

This work introduces a similarity solution to the problem of a viscous, incompressible and rotational fluid in a right-cylindrical chamber with uniformly porous walls and a non-circular cross-section. The attendant idealization may be used to model the non-reactive internal flow field of a solid rocket motor with a star-shaped grain configuration. By mapping the radial domain to a circular pipe flow, the Navier–Stokes equations are converted to a fourth-order differential equation that is reminiscent of Berman’s classic expression. Then assuming a small radial deviation from a fixed chamber radius, asymptotic expansions of the three-component velocity and pressure fields are systematically pursued to the second order in the radial deviation amplitude. This enables us to derive a set of ordinary differential relations that can be readily solved for the mean flow variables. In the process of characterizing the ensuing flow motion, the axial, radial and tangential velocities are compared and shown to agree favourably with the simulation results of a finite-volume Navier–Stokes solver at different cross-flow Reynolds numbers, deviation amplitudes and circular wavenumbers.

Self-similarity solution of NS equations for porous non-circular tube 399 2. Solution methodology 2.1. Geometry and governing equations Our geometric configuration consists of a semi-infinite porous cylinder that is closed at the headwall while admitting a uniformly distributed sidewall injection. Considering any real function α with zero mean, the cross-section is determined by the tangentially evolving radius, r w = a(1 + εα), where a is fixed and ε is a parameter. As shown in figure 1, the fluid is injected perpendicularly to the wavy surface at a constant speed U w . In view of this geometry, the steady incompressible Navier-Stokes equations are considered with all spatial coordinates, time, velocity and pressure variables made dimensionless using a, a/U w , U w and ρU 2 w , respectively. The corresponding cross-flow Reynolds number may be written as Re = U w a/ν, where ν denotes the kinematic viscosity. Then using u = (ur, uθ , uz) and p for the velocity and pressure, our system of equations reduces tō ∇ · u = 0 and u ·∇u = −∇ p + 1 Re∇ 2 u, (2.1a,b) where the overbar in∇ denotes differentiation with respect to the unmodified dimensionless variables (r,θ,z).

Mapping to a circular cross-section
To simplify the application of the boundary conditions at the sidewall, and thus avoid the need for Taylor-series expansions about a fixed radius, it is helpful to map the dimensionless spatial coordinates (r, θ, z) in such a manner to relocate the wall boundary conditions to a spatially invariant point. This may be accomplished through the coordinate transformation: , θ = θ , z = z. (2. 2) The vector basis (e r , e θ , e z ) defined in figure 1 remains unchanged. Although this substitution affects the form of the radial and tangential derivatives, the wall surface can now be specified using A(r) = r − 1 = 0. The corresponding normal unit vector, n = n r e r + n θ e θ , becomes and so the wall-normal injection condition simplifies to u = −n at r = 1. Naturally, this substitution gives rise to a more complex form of the Navier-Stokes equations, which is detailed in appendix A.

On the minimum perturbation order
In their elegant analysis of the axisymmetric flow in a porous tube, Balachandar et al. (2001) establish a framework in which small perturbations of the form ε cos(mθ ) may be imposed at the sidewall. By examining the asymptotic behaviour of the ensuing solution as r → 0, a viscous core is shown to be present near the centreline. Accordingly, viscous effects dominate on the scale of r = O(1/ √ Re) for wall perturbations of O(ε). Furthermore, as we enter the viscous core region, the 400 M. Bouyges, F. Chedevergne, G. Casalis and J. Majdalani a FIGURE 1. Sketch of the geometric configuration and coordinate systems corresponding to a non-circular, wavy cross-section with α = cos(mθ). The sketch corresponds to m = 7 lobes and a 20 % radial deviation in ε.
axial vorticity, which accompanies the evolution of a finite tangential velocity, appears at O(εRe). Consequently, the linear perturbative analysis remains valid as long as εRe 1. To overcome this practical limitation, Balachandar et al. (2001) introduce a nonlinear patch on the scale of r = O( √ ε). This viscous correction leads to a solution that extends the range of applicability to εRe 1, where the axial vorticity increases to O(1). A similar conclusion may be reached in our analysis, which requires second-order perturbations of O(ε 2 ) to adequately capture the nonlinear interactions that evolve at previous orders. Moreover, a second-order approximation will be warranted to satisfy mass conservation requirements by virtue of the angular dependence being of the harmonic form, e imθ . The approach that we follow is therefore compatible with Balachandar's framework, as it seeks to derive a second-order viscous solution for a periodically distorted cross-section. In this process, however, an effort will be made to secure the no-slip condition in all three spatial directions, including the axial and tangential velocities at the sidewall, thus leading to a uniformly valid solution that extends from r = 0 to 1 inclusively.
To justify the need to carry out a second-order approximation, it may be instructive to revisit the mass flow rate evaluation at the sidewall. We start by remarking that, for an identical mean radius, the circumference of the cross-section corresponding to a wavy wall, such as a star-shaped grain configuration, will always exceed that of its counterpart with a fixed circular radius. Then, recognizing that α is a periodic function with a vanishing mean value, the mass surplus can be estimated using an asymptotic expansion ofṁ w (z), namely, The resulting equality confirms that the first-order corrections have no bearing on the injected mass flow rate. In particular, the axial velocity, which contributes to the mass injection, cannot by itself secure the mass balance without taking into account its interactions at the second order. Moreover, a second-order representation will enable us to extend the validity of the model to larger deviation amplitudes with ε 2 1, while simultaneously aiding in improving the permissible range of Reynolds numbers to εRe 1, as predicted by Balachandar et al. (2001).
Self-similarity solution of NS equations for porous non-circular tube 401 2.4. Asymptotic expansion Given the above considerations, the velocity and pressure fields, (u, p), may be decomposed using where superscripts denote successive asymptotic orders. This two-term approximation may be substituted back into the Navier-Stokes equations, expressed through (A 1), and then further linearized using binomial expansions of (1 + εα) −1 terms. These operations enable us to identify and segregate three sets of equations that appear at O(1), O(ε) and O(ε 2 ), consecutively. A similar expansion may be applied to the wall injection boundary condition (2.3), thus leading to When translated to the asymptotic velocity corrections, these expressions give rise to an assortment of conditions that must be imposed at successively increasing orders, namely, The solution to the three sets of equations and their associated boundary conditions will be described in § § 2.5-2.7. Furthermore, to justify the use of the upcoming similarity transformation, it will be helpful to examine the mass conservation requirement. By equating the radial mass flow rateṁ w (z), which originates from the sidewall over a length of tubing that extends from the headwall to a point z, to the axial mass flow rateṁ cross (z), crossing section z, one arrives at the following equality: A sufficient condition for u z to satisfy the previous relation (2.8) is to assume a linear variation with respect to z, as foretold by Hiemenz (1911). To make further headway, only separable product solutions are sought, thus turning the linear dependence on z into a necessary condition. Procedurally, this enables us to specify the axial velocity corrections, u (0) z , u (1) z and u (2) z , with a linear dependence on z.
2.5. Leading-order analysis In the absence of angular wall deformation (ε = 0), the classical problem of an injection-driven motion in a circular cylinder is restored. Being independent of θ , one recovers u (0) θ = 0 and ∂ θ f = 0 for any scalar function f . In this situation, Berman's formulation may be employed with a streamfunction ψ(r, z) = zF(r), where F(r) 402 M. Bouyges, F. Chedevergne, G. Casalis and J. Majdalani represents the characteristic mean flow function (Berman 1953). The leading-order solution may be subsequently retrieved from Berman's fourth-order nonlinear ordinary differential equation: In the above, the boundary conditions may be ascribed to axisymmetry and wall-normal injection. In hindsight, the solutions obtained by Taylor (1956) and Culick (1966) correspond to Berman's formulation in the limit of Re → ∞. In fact, the Taylor-Culick profile begins to resemble Berman's solution for Re > 100. Nonetheless, as we aim to extend the validity of the viscous model to a generic non-axisymmetric configuration, we continue to employ Berman's model and retain its dependence on the cross-flow Reynolds number. Thus, for a given Re, (2.9) may be solved numerically for F(r), and so the three velocity components and pressure at order O(1) may be deduced directly from where both P 0 and κ 0 represent undetermined integration constants that may be obtained from Berman's expression in the modified coordinate system (r, θ , z).
2.6. First-order analysis The system of linear equations (A 1) is written at the order O(ε) and rearranged into (C 1). Guided by the wall-normal injection conditions requiring the appearance of α and its derivative, a solution is sought using separation of variables. Note that no additional axial pressure gradient is required for product solutions, since the mass flow surplus (2.4) is of second order. The solution then simply reads The substitution of these expressions into (C 1) gives rise to a θ-dependent system in which every equation can be collapsed into the form of φ(r) = −α (θ )/α(θ ), where φ denotes an operator that may be constructed from the aforementioned functions and their derivatives. Naturally, without any loss of generality in the space of real functions, harmonic variations of the form α(θ ) = cos(mθ ) are retained, m being a constant. Since m represents the number of lobes along a fixed circular circumference, it will be referred to hereafter as the circular wavenumber. The underlying consequence of this choice for α is that the mass flow conservation can be expressed at the third order usinġ For the reader's convenience, the θ-dependent system of equations is provided in (C 2). The system consists of one first-order and three second-order ordinary differential equations. The corresponding injection conditions, which stem from (2.7), lead to v θ (1) = 1 and v r (1) = v z (1) = 0. In order to achieve closure, four additional constraints are still necessary, and these may be specified at r = 0 by introducing Taylor-series expansions of the different functions with respect to r in (C 3). After some effort, we arrive at v r (0) = −F (0)/2 and v r (0) = v θ (0) = v z (0) = 0, with the additional realization that m 2 = (0, 1), to avoid a mathematically incongruent outcome. In practice, we are only interested in cross-sections with m 3.
With the advent of a sufficient number of boundary conditions, a collocation method, such as the one described by Canuto et al. (1988), may be used to discretize the system of equations and recover the first-order velocity and pressure corrections.

Second-order analysis
To make further headway, we proceed to decompose the second-order corrections using the following classical forms: On the right-hand side of the resulting momentum equations at O(ε 2 ), terms that depend on the two previous orders, which are caused by nonlinear interactions, may be explicitly seen (appendix C). Computing this second-order correction will therefore enable us to capture the nonlinearities alluded to by Balachandar et al. (2001) when εRe is no longer a small quantity. The resulting set of seven linear equations can be subsequently consolidated into two concise and independent systems: the first assortment consists of four linear equations that prescribe (u (2) 2 , p (2) 2 ), as given by (D 1), and the second entails three linear equations that control (u (2) 0 , p (2) 0 ), as given by (D 2). Being independent, these two systems may be solved separately. Moreover, since u (2) θ only appears in the first system, each set leads to a well-posed problem with the same number of equations and unknowns, except for the constant κ. Its determination will be described separately.
The boundary conditions associated with u (2) 2 may be readily obtained from (2.7), specifically at the wall and the centreline, using the r → 0 Taylor-series expansions introduced previously. These produce (2.14) With the boundary conditions in hand, the computation of (u (2) 2 , p (2) 2 ) may be achieved easily using a discrete collocation, for instance.
Unlike (u (2) 2 , p (2) 2 ), the three linear equations for (u (2) 0 , p (2) 0 ) may be carefully reduced into a single third-order ordinary differential equation that depends solely on u (2) r,0 . This equation is given by (D 3). The undetermined pressure constant κ, which appears on the right-hand side of (D 3), has a direct influence on u (2) z,0 (1). Since κ appears as a forcing term in a linear relation, the connection between u (2) z,0 (1) and κ proves to be linear. At the outset, a unique value for κ may be retrieved for a given set of input parameters in order to secure u (2) z,0 (1) = 0. For the reader's convenience, its values for different ranges of Re and m are provided in table 1.

Characterization at successive orders and wavenumbers
In summary, the velocity and pressure fields may be expressed using (2.15) To illustrate the behaviour of the solution, we confine our attention to a polar slice that extends between θ = 0 and 2π/m. In figure 2, we examine the sensitivity of the solution for u r and u z at increasing asymptotic orders using characteristic values of ε = 0.1, Re = 100 and m = 7. Only the radial variation is captured here, which may be accomplished by plotting the solution along a spoke at θ = 0 and a fixed value of z. The radial velocity, which is featured on the left-hand side of the graph at three increasing orders of ε, confirms that even the first-order approximation, u (1) r , is considerably dissimilar from the leading-order solution. This disparity implies that the radial displacement of the sidewall leads to a significant distortion of the mean flow field compared to the circular configuration. However, these results also suggest that the second-order correction, u (2) r , only weakly affects the overall radial velocity. Although similar observations may be inferred for the tangential velocity, they are omitted here for the sake of brevity.
In contrast, by turning our attention to the axial velocity variation on the right-hand side, it may be seen that the second-order correction leads to a non-negligible shift in u z , especially as the centreline is approached, where the first-order contribution is nil by virtue of u (1) z (0) = 0. The need for a second-order correction is therefore pivotal. In fact, the discrepancy near the centreline between Berman's leading-order solution and that at O(ε 2 ) may be attributed to the mass flow rate increase reported in (2.4). Recalling that the mass flow rate entails an integration over θ, the contribution of the first-order correction u (1) z multiplied by the periodic function cos(mθ ) vanishes identically. However, the contribution of the second-order correction u (2) z , which Self-similarity solution of NS equations for porous non-circular tube involves a θ-independent term (u (2) z,0 ), leads to a finite contribution. Based on this simple comparison and mathematical verification, it may be concluded that the first-order correction affects all three velocity components, thus leading to marked differences from the circular configuration. As for the second-order correction, it proves to be essential to properly capture the behaviour of the axial velocity, especially near the centreline, and therefore helps to secure the principle of mass conservation.
To better understand the effect of m on the second-order approximation, both u θ and the centreline value of u z are depicted in figure 3 at ε = 0.1 and Re = 100. The tangential velocity is featured in figure 3(a) along a spoke with θ = π/(4m) and m = 5, 6, 7, 8 and 9. Interestingly, the effect of the injection condition in (2.7) on the magnitude of u θ may be captured at the sidewall, where |u θ | increases incrementally with m. The maximum values of u θ also shift outwards towards the sidewall with successive increases in m.
It is also useful to explore the sensitivity of the centreline value of u (2) z (0) to variations in m ∈ [3, 12] in figure 3(b). Since u (2) z,2 (0) = 0, it is sufficient to illustrate the behaviour of u (2) z,0 (0) in order to infer the underlying trend. Although the relation between u (2) z,0 (0) and m remains concealed in (D 2), the corresponding curve suggests the existence of a monotonically increasing quadratic polynomial of the form u (2) z,0 (0) = c 0 + c 1 m + c 2 m 2 , whose coefficients in figure 3(b) are c 0 = −1.242, c 1 = −0.1875 and c 2 = 0.817. Despite the appearance of m 2 in the equations explicitly, determining the quadratic solution analytically for any Re and ε proves to be a daunting task. Nonetheless, based on this characteristic graph, it may be ascertained that the role of u (2) z,0 (0) can be significant, especially at increasing values of m. Although not shown, we also find that the effect of increasing the Reynolds number on u (2) z,0 (0) is negligible in comparison to m. As for flow rotationality, figure 4 displays isocontours of the tangential component of the steady-state vorticity Ω z in an arbitrary (r, z) plane; also provided is a comparison of the streamlines for the present solution side by side with Berman's. Based on this graph, it may be seen that the star-shaped configuration leads to deeper penetrating streamlines than the circular case. Specifically, the solid lines, which correspond to the present formulation, consistently approach the chamber axis more rapidly than the dashed lines of the non-wavy solution.

Comparison to Navier-Stokes simulations
To further illustrate the behaviour of the non-circular profile and more accurately define its domain of applicability, a comparison is now undertaken between the results obtained from the present semi-analytical formulation and those predicted by a Navier-Stokes solver. To this end, the Navier-Stokes solver is introduced first along with its computational characteristics. Second, the sensitivity of the flow computations to variations in the three main parameters, ε, Re and m, are explored and compared to our second-order approximations.
3.1. Solver description Our numerical simulations are performed using a compressible Navier-Stokes solver called CHARME, which belongs to the multi-solver suite CEDRE developed by the French Aerospace Laboratory (ONERA). A detailed overview of the code's functionality is provided by Refloch et al. (2011), while examples of its computational capabilities, which include laminar flow simulations and verifications of linear stability analyses of rocket internal flow fields, are furnished by Chedevergne et al. (2012) Self-similarity solution of NS equations for porous non-circular tube 407 and Boyer et al. (2013). In these studies, an excellent agreement is reported between theory and computations, thus showcasing the solver's effectiveness in capturing the mean flow development as well as the unsteady hydrodynamic instabilities that evolve in a rocket chamber. As usual, the latter may be idealized as a porous tube with uniform wall injection, which is consistent with the geometric configuration used presently. As for the code, a second-order discretization stencil is implemented for spatial interpolation, while Euler's fluxes are computed using a Roe scheme. Throughout these simulations, the length and mean radius of the porous tube are kept constant at 8a and a = 1 cm, respectively. Although the value of a may be modified, the aspect ratio is held constant. Furthermore, the wall injection speed U w is adjusted to calculate the targeted Reynolds numbers while ensuring that it remains sufficiently small to promote low compressibility levels. At the headwall, a slip condition that is compatible with our model is applied to avoid the development of an undesirable boundary layer. In the domain exit plane at z = 8a, a pressure outlet condition is imposed. The code is fully three-dimensional and the assumption of linear axial variation is not forced. Then, using a three-dimensional unstructured mesh with approximately four million total elements, 600 cells are selected to construct the axial grid along the sidewall, whereas 150 cells are distributed radially along the diameter. Although the mesh size can generally depend on Re, ε and m, the present discretization size is confidently chosen because of its demonstrated effectiveness in previous studies of analogous flow fields (Chedevergne et al. 2012;Boyer et al. 2013). In what follows, the term 'model' will be used to denote second-order asymptotic results obtained from the present framework, whereas 'CFD' will refer to numerical (computational fluid dynamics) simulations acquired from the nonlinear Navier-Stokes solver.

3.2.
On the linearity of the axial velocity with respect to z Inspired by Berman's solution and the mass conservation relation (2.8), the axial velocity component can be assumed to be linear with respect to the axial coordinate z. In order to verify that the CFD solution will indeed mimic this linear behaviour, figure 5 is used to compare our model to the computed solution for ε = 0.1, Re = 100 and m = 7. To magnify the possible departures from the expected linearity, the leading-order contribution of O(1), i.e. Berman's solution, is subtracted from the axial velocity component u z , thus leaving only first-and second-order corrections. Based on this comparison, it is clear that the computed axial velocity remains linear with respect to z up to approximately z = 7 using r = 0.3. Beyond this location, it is possible for the outlet pressure condition to slightly influence the axial development of the flow field. Apart from this outlet effect, it may be argued that the Navier-Stokes solver in figure 5 helps to validate the assumption of a linear variation with respect to z. Moreover, the visual agreement in the resulting slopes confirms that the simulation accurately reproduces the radial and tangential evolutions of the modelled flow. In fact, additional comparisons only reinforce this behaviour. In the interest of simplicity, subsequent results will be displayed halfway through the chamber at z = 4.

Flow characterization based on both asymptotics and computations
To take advantage of the periodicity in the angular direction, flow field comparisons are reproduced in a polar slice that extends from θ = 0 to 2π/m. Note that the mesh used in conjunction with the Navier-Stokes simulations does not rely on periodic conditions; the tangential coordinate θ extends over the full range of [0, 2π]. We begin 408 M. Bouyges, F. Chedevergne, G. Casalis and J. Majdalani in figure 6 with the isocontours of the three velocity components using ε = 0.1, Re = 100, z = 4 and m = 7. Despite the linearization that affects the second-order asymptotic framework, an excellent agreement may be noted between the isocontours generated from the numerical simulations and those obtained from the asymptotic model for all three radial, tangential and axial velocity components.
By turning our attention first to the tangential velocity, it may be confirmed that its maximum absolute values occur at θ = π/(2m) and 3π/(2m) taken along the sidewall because of surface waviness; as for its sign, it depends on whether the flow is situated above or below the meridian line at θ = π/m, where the tangential velocity vanishes identically. The tangential velocity also vanishes at the domain delimiting lower and upper spoke lines located at θ = 0 and 2π/m, where the radial velocity u r matches the wall-normal velocity in both direction and magnitude. At those specific sites, both axial and tangential velocities vanish, thus leaving the flow to enter the chamber radially inwards towards the centreline (i.e. with no directional distortion). Naturally, |u θ | diminishes while moving inwards away from the sidewall, or when moving tangentially towards the meridian lines as well as the upper and lower domain delimiting lines. In anticlockwise fashion, u θ switches from negative to positive values as we cross from 0 < θ < π/m to π/m < θ < 2π/m. Moreover, given that the tangential velocity is induced by small deviations from a circular cross-section, it remains smaller than its radial and axial counterparts, being asymptotically driven by O(ε) surface displacements, unlike the radial and axial velocities, which are chiefly induced by the mass injection mechanism.
At this juncture, it may be helpful to recall that, in the absence of surface waviness, the isocontours of the axisymmetric axial velocity lead to concentric circles that exhibit a maximum value along the centreline. In the present configuration, the isocontours of u z stretch radially outwards along the meridian line of θ = π/m, where the chamber radius happens to be the shortest. Naturally, the area contraction in flow cross-section at this angular position induces a local increase in the axial velocity. Conversely, these isocontours shift radially inwards at the upper and lower boundaries where local wall expansions lead to an increase in the flow cross-section and, therefore, a corresponding reduction in u z .
Self-similarity solution of NS equations for porous non-circular tube In fact, continuity may be used to partly explain the trends observed in the isocontours of the radial velocity. Upon close inspection of the corresponding graphs, it may be ascertained that u r vanishes along the centreline, where u z reaches its peak value, and that it increases in regions of axial velocity deficit. This behaviour may be attributed to the flow turning mechanism that affects the wall-injected stream, specifically as it negotiates a 90 • turn before merging axially with the core flow. At θ = π/m, owing to the stretched axial velocity isocontours, the radial velocity deceleration towards the centreline is more rapid than at the θ = 0 and 2π/m borderlines, where the flow may be seen to be radially expanded. Another characteristic of u r , which has been well documented in the literature, corresponds to its velocity overshoot compared to its wall value (Saad & Majdalani 2010). Despite the unitary speed of |u r | = 1 at the sidewall, |u r | undergoes first an increase in absolute value due to the sudden radial pinching of the circumferential area that is normal to the injected flow, and which does not permit the flow to develop a sufficient u z component to adequately transport the mass axially. The resulting increase in |u r | may be noted in the peak contour points in figure 6, namely, those that materialize at some intermediate positions between the wall and the centreline. Naturally, when the flow is pinched, the locus of these extrema moves closer to the centreline. As alluded to earlier, another byproduct of waviness and viscosity may be the production of axial vorticity, which is absent in the strictly axisymmetric Taylor-Culick model (Balachandar et al. 2001). To illustrate the corresponding behaviour, isocontours of axial vorticity, namely, Ω z , are displayed in figure 7 for ε = 0.1, Re = 100, z = 4 and m = 7. Using the variable coordinate transformation at the wall, Ω z may be expressed as (3.1) Note that this vorticity component vanishes for the axisymmetric case. In order to improve the velocity gradients and thus the vorticity, the spatial integration scheme within the Navier-Stokes solver is increased to the third order in figure 7. Graphically, one may infer that the Navier-Stokes computations tend to agree qualitatively with the second-order asymptotic solution, especially in the two bulk regions where the axial vorticity is largest. Nonetheless, discrepancies may also be seen at the sidewall near θ = π/(2m) and 3π/(2m). The bulk vorticity regions that develop in this configuration are consistent with those detailed by Kurdyumov (2006). Furthermore, the dissimilarities observed near the sidewall may be attributed to the practical limitations associated with asymptotic expansions. They may also be ascribed to the compounded effects of numerical errors that accrue during the evaluation of velocity gradients within the Navier-Stokes solver.
3.4. Sensitivity to various Reynolds numbers, radial deviations and wavenumbers To further explore the effects of varying the Reynolds number, radial deviation amplitude and circular wavenumber on the flow development, the sensitivity of the three velocity components to these parameters may be quantified by comparing the second-order corrections to the individual computations taken along different polar angles. Guided by the isocontours of figure 6, both u r and u z seem to display interesting structures along angular cuts taken at θ = 0 and π/m, whereas u θ may be best featured at θ = π/(2m). We therefore proceed to evaluate the three velocity  (2) z ]/z and u (1) r + εu (2) r with m = 7 crests.
components along their most representative inclination angles and showcase their outcomes in figures 8-10 using both asymptotics and computations.
As we seek to better understand the role of ε on the axial velocity, the second-order corrections obtained asymptotically are compared to those computed in figure 8(a) using Re = 100, z = 4, m = 7 and two values of ε that correspond to 0.01 and 0.1. Results are deliberately displayed along a constant θ = π/m spoke where the largest deviations from the axisymmetric motion may be realized. To further magnify differences, the exact Berman solution that is recovered at leading order is subtracted from both solutions, and the outcome is divided by ε, thus leaving us with [u (1) z + εu (2) z ]/z. In view of the dual magnification levels that are implemented, it may be argued that the agreement between computations and asymptotics is satisfactory, especially in the ability of the model to capture both sidewall and centreline variations with ε rather consistently. The most noticeable discrepancies remain small and mostly confined to the annular region that extends approximately between r ≈ 0.4 and 0.7.
In fact, because of the magnification effect, these discrepancies actually fall within the numerical uncertainty that accompanies our simulations. The reasons are these. Since the mesh is uniformly discretized along the z direction, the accuracy of u z /z directly depends on the finite gap between two axial positions. Furthermore, the disparities can stem either from higher-order terms that are not accounted for in the asymptotic formulation, or from the diffusive nature of the second-order discretization scheme that is adopted in the Navier-Stokes solver.
As for the influence of the Reynolds number on the radial velocity, it is illustrated in figure 8(b) using ε = 0.01, m = 7 and two values of Re that correspond to 100 and 1000, respectively. Choosing a relatively small value of ε enables us to compute the flow field with a large Reynolds number. Here, too, we display only the secondorder correction in the radial velocity, u (1) r + εu (2) r , by subtracting the leading-order contribution. The agreement that is achieved between the numerical simulations and the asymptotic formulation is gratifying to note, especially as the Reynolds number is set at 1000.  FIGURE 9. Sensitivity of the tangential and axial velocities to variations in m using both computations (symbols) and second-order asymptotics; here we take ε = 0.1 and Re = 100. Results are shown along constant angular spokes of (a) θ = π/(2m) for the tangential velocity u θ /ε = u (1) θ + εu (2) θ and (b) θ = 0 for the axial velocity [u (1) z + εu (2) z ]/z.
In figure 9(a), the sensitivity of the tangential velocity on the circular wavenumber is characterized by showcasing u θ /ε = u (1) θ + εu (2) θ along a polar angle of θ = π/(2m) using ε = 0.1, Re = 100 and m = 5, 7 and 9. The skipping of even values of m reduces visual clutter by helping to distinguish between the remaining curves, especially since the use of m = 4, 6 and 8 leads to nearly identical trends. It is also guided by practical reasons, such as the preferred grain perforations in industry, which are invariably designed with an odd number of star points, crests or lobes, with the hope that such arrangements would help to mitigate the development of acoustic instabilities. Here, too, the agreement between computations and the modelled solution appears to be excellent in the core region. However, small differences begin to appear when the wavenumber is increased or when the sidewall is approached. We also remark that at θ = π/(2m), the contribution of second-order terms to the tangential velocity become negligible because of the sin(2mθ) term that appears in the definition of u (2) θ . As one may surmise from (2.13), the minor disparities that are detected near the sidewall may be attributed to the omission of third-order terms in the asymptotic expansion of the injection velocity. Asymptotically, one must have where the quantity α (α 2 − 1 2 α 2 )ε 3 reduces to m 3 ε 3 /2 at θ = π/(2m); this third-order correction actually matches the gap between the computational and second-order asymptotic values at the wall, which widens with successive increases in m and ε.
The sensitivity of the axial velocity on the circular wavenumber is illustrated in figure 9(b) along a different polar angle, namely θ = 0, at the same Reynolds number and deviation values. Here, Berman's contribution u (0) z is subtracted from the total axial velocity component in order to showcase the correction only, as done previously in figure 8. Consistently with the behaviour observed in the tangential velocity comparison, the agreement between computations and the modelled solution appears to be quite satisfactory for m = 5, although it slowly deteriorates as m is increased. As before, the overshoot of the model near r ≈ 0.9 may be attributed to the diffusivity of the Navier-Stokes solver; the accuracy of the simulation along the axis remains limited by the finite discretization in the axial direction. Since the absolute value of the axial velocity along the axis |u (2) z (0)| increases with m, so does the uncertainty associated with the axial mesh discretization, thus leading to a growing disparity between the simulations and the modelled solution. Furthermore, it may be helpful to note that the displayed velocity components are divided by ε to better differentiate among the curves. Nonetheless, the noted discrepancies remain small compared to the total velocity magnitude u.
Before leaving this discussion, it may be instructive to note that, in the foregoing work, the product εRe has been consistently smaller than 10. Nonetheless, comparisons could have been undertaken for values that reach εRe ≈ 30. Above this level, the Navier-Stokes solver will have difficulty converging onto a steady periodic solution. An example of the model's validity for εRe = 20 is provided in figure 10, where isocontours of the velocity are displayed at two different values of m. By comparing the upper and lower sections of the graph, the ability of the model to mirror the computed solution seems gratifying, except in what concerns the negative peak magnitudes of u r in the model, which overshoot their computed values. The axial velocity in figure 10(d) also exhibits different lobe shapes when comparing the modelled solution to its simulation at the same contour levels.

Conclusion
In this study, the Navier-Stokes equations are applied to the asymmetric problem that arises in the context of an injection-driven mean flow motion in a porous tube with a wavy non-circular cross-section. After mapping the radial coordinate to a fixed 414 M. Bouyges, F. Chedevergne, G. Casalis and J. Majdalani boundary, the equations of motion are transformed and shown to recover Berman's fourth-order differential equation in the absence of surface distortion. To make further headway, an asymptotic expansion is pursued in the deviation amplitude, ε, which stands for the maximum radial displacement in a unit circle. This enables us to derive the system of equations that prescribe the stationary motion both at the first and second orders in ε. Then using judiciously imposed similarity transformations, the first-order system is reduced to four ordinary differential equations that can be solved straightforwardly. Along similar lines, the second-order system is transformed into an assortment of seven ordinary differential equations that can be solved with equal ease. The second-order formulation is subsequently shown to provide a practical similarity-based, asymptotic approximation to the problem at hand. Comparisons to nonlinear numerical simulations serve to verify the accuracy of the solution and to confirm the necessity of carrying out the analysis at least to second order. In this vein, we find that the axial velocity component is strongly affected by the second-order correction, which also plays a key role in securing mass conservation, especially in long chambers. From a fundamental standpoint, a high-order approximation not only extends our range of applicability to larger values of the Reynolds number and deviation amplitudes, it also helps to capture the axial vorticity that evolves in the chamber, as well as the viscous patch that develops along the centreline. Moreover, the semi-analytical solution enables us to assess the sensitivity of the flow field to variations in the cross-flow Reynolds number, the deviation amplitude and the circular wavenumber, which is specified by the number of lobes or crests in an actual grain perforation.
The inception of a viscous second-order mean flow model as a substitute to the Taylor-Culick profile opens up a parallel line of research inquiry into the stability and performance of rocket chambers in the presence of wall distortion. Although numerous stability investigations of the Taylor-Culick profile exist, with and without particle mean flow interactions, the corresponding problems in non-circular grain configurations remain relatively unexplored. In a circular motor, for example, the role of the tangential vorticity component on stability is relatively well understood. However, in the presence of asymmetry and surface distortions, the impact of axial vorticity on instability remains an open question. Given the deeply modified vorticity structures in a non-circular configuration, it will be interesting to explore the effects of deviation amplitudes and wavenumbers on flow instability. It is also hoped that these and other research questions will be addressed in future work.
Thus the system of equations that controls v r , v θ , v z and q may be rearranged into: with the auxiliary boundary conditions v z (1) = v r (1) = 0, v θ (1) = 1, (C 3a,b)