The Role of Buoyancy Induced Instability in Transpirational Cooling Applications

Transpirational cooling is an effective thermal protection method in hypersonic vehicles. In order to properly manage the high heat load, an understanding of the convective flow regimes as well as the thermophysical properties of the working fluid are required. Often, the vehicle’s fuel is re-purposed as the coolant or working fluid that is passed through the porous media. If the geometry is such that the coolant is heated from below, buoyancy-induced instability can ensue resulting in a mixed convection phenomena. Transpirational cooling applications require a unique analysis which combines a Darcy–Forchheimer relationship for the momentum relation, a flowing base state which introduces non-negligible convective terms for the energy equation, and a novel consideration of a cubic density dependence on temperature. This latter feature is justified by fitting thermodynamic data for typical transpirational cooling conditions. A base state solution is provided and the onset of instability is investigated using linear stability analysis. The governing equations are solved utilizing multiple methods, comparing results from a combination of analytical solutions, finite difference, power series, and Chebyshev methods. Results demonstrate excellent consistency in predictions across these methods and indicate that including non-linear density effects promote a stabilizing effect. Finally, the effect of varying the net through-flow in the porous media is investigated.


Introduction
Transpirational cooling is an effective thermal protection method in hypersonic vehicles. The transpirational cooling concept was first applied to the throat of rocket engines in the 1940s. Since then it has been applied to other parts of the aircraft that experience high heat loads such as the leading edge and fuel injection strut. Often, the vehicle's fuel is re-purposed as the coolant or working fluid that is passed through a porous media section to form a continuous thin film over the protected wall. Thus, a proper understanding of the underlying mechanisms in transpirational cooling applications leads to more robust design capable of handling higher heat loads.
Early investigation of transiprational cooling phenomena was carried out by Rannie [1] who analytically modeled the porous wall temperature to investigate the laminar transpired sublayer. His model aligned with experimental data at low thermal conductivity for the porous substrates. Other analytical approaches for transpirational cooling have also been detailed by Eckert [2,3] including how to model coolant injection and hot gas heat transfer from an exterior hot boundary layer. Recently, a closed form analytical solution was developed for 2D pressure and velocity field for transpirational cooling applications [4]. Several experimental works have been put forth detailing the effects of the transpirant, mass flow rate, inlet temperature, flow coking, etc. on cooling efficacy and performance [5][6][7][8]. Dahmen et al. [9] analyzed the experimental transpirational cooling setup in Ref. [6], using numerical simulation to provide full flow field details that were difficult to acquire through experiments. Other computational studies have been performed that quantify cooling efficacy versus inlet conditions, thermal properties, and boundary conditions [10,11].
All previous studies have either analytically, experimentally, or numerically, investigated transiprational cooling in geometric configurations where the hot surface is on the top or side and the cold surface is on the bottom or side. However, if the geometry is such that the hot surface and cold surface are on the bottom and top, respectively, then Bénard-type convection can ensue resulting in a mixed convection phenomena. This type of arrangement occurs on the bottom wall of the aircraft or the top wall of the combustion chamber. In these cases, an understanding of the buoyancy-induced instability mechanisms need to be described and understood. The tools of linear stability analysis provide a pathway for this understanding.
A similar geometric and flow configuration to that occurring in transpirational cooling applications, albeit with no base flow, is the Horton-Rogers-Lapwood (HRL) problem. The HRL problem consists of a buoyancy-induced instability in porous media due to heating from below [12,13]. The HRL problem is extensively described and characterized in Nield and Bejan [14]. A pioneering variation of the HRL problem where base flow is perpendicular to the base temperature gradient was studied by Prats [15]. Several variations of the HRL and Prats problems have been investigated by . The variation of the HRL problem most relevant to the transpirational cooling problem at hand is with a base flow parallel in the direction of the base temperature gradient, i.e., vertical base flow. This problem was first studied by Wooding [19] on a semi-infinite (vertical) porous domain. Sutton then characterized stability in a finite domain-of both width and height-against through-flow velocities of a small magnitude [20]. Initially, a small negative through-flow has a slight destabilizing effect which transitions to a stabilizing effect as velocity is increased. This effect is exaggerated as the domain width becomes smaller. Next, Homsy and Sherwood investigated the stability theory on a finite domain for both strong and weak through-flows [21]. They found that mass discharge in the positive vertical direction destabilizes while the opposite was true for negative through-flow.
The previous works provide much insight to buoyancy-induced instabilities occurring in porous media applications such as transpirational cooling. However, transpirational cooling has additional complexities which are not accounted for by the hydrodynamic stability analysis in these works. These include the highly non-linear density variation with temperature, i.e., a nonlinear Boussinesq approximation [22,23]. Due to the often extreme temperature difference between the inlet and exit of the porous media, extremely large density gradients can occur. These density gradients arise from heating cryogenic fuels like methane or hydrogen to a supercritical state. Real-gas compressibility effects can be significant at these extreme temperatures and pressures. Thus, the objective of this work is to include the effects of a non-linear cubic density dependence on temperature in a base state configuration with a vertical through-flow to quantify and understand the instability mechanisms occurring in the porous medium. Three cases are considered in this work, each with an additional layer of complexity: (1) A traditional HRL configuration, (2) a HRL configuration with a non-linear density variation in temperature, and (3) adding a base vertical through-flow to Case 2. The third case is most representative of transpirational cooling. The non-dimensionalization introduced in this study to account for the drag term results in a Grashof-type scaling rather than the typical Rayleigh-type scaling encounter in HRL (or Wooding) style problems. The resulting equations are solved using a power series method, Chebyshev method, and finite difference method. The three cases of increasing complexity elucidate how the various effects of a non-linear density and a varying vertical through-flow affect the stability in transpirational cooling systems.
This work is divided into five sections. First, the mathematical framework used to study transpirational cooling is described. Second, linear stability analysis is performed and detailed. Next, an overview of the various solution techniques are outlined. Finally, results are given followed by a discussion and conclusions.

Governing Equations
A simplified model for transpirational cooling is presented in Figure 1. Here, a fluidsaturated isotropic porous medium of infinite horizontal extent is confined between two permeable boundaries a distance d apart. The hot gas flowing over bottom layer is modeled as a constant temperature boundary condition as in Ref. [6]. Coolant continuously flows over the top permeable boundary (z = d in Figure 1), and creates an isothermal condition. Wall suction provides a constant coolant injection aligned with the gravitation vector [3]. This following problem configuration is similar to the Wooding problem [19], where a constant through-flow is heated which gives rise to an additional thermally-induced natural convection of the Beńard type. These assumptions, along with the Boussinesq approximation, lead to Dirichlet conditions for velocity and temperature at the top and bottom boundary. Analysis of this problem requires conservation of the mass, momentum, and energy equations subject to the aforementioned boundary conditions. In regards to the momentum equation, several experimental [6,7] and numerical studies [8,9] have determined that transipirational cooling applications are accurately approximated by the Darcy-Forchheimer equation. The curl of the Darcy-Forchheimer equation is presented here to drop the pressure term. The energy equations for the porous material and the fluid phase can be simplified to one equation if local thermal equilibrium exists between the solid and fluid phases, and there is no net conduction between the phases [14]. The resulting single equation uses an overall heat capacity and thermal conductivity per unit volume of the porous material and fluid denoted by the subscript 'm'. For example, (ρc p ) m = (1 − ψ)(ρc p ) s + ψ(ρc p ) f , where ψ is the porosity and the subscripts 's' and 'f' refer to the solid and fluid phase, respectively. Langener and Wolfersdorf [6] investigated these assumptions along with negligible work due to pressure, negligible viscous heating, and constant density in the energy equation to study transpirational cooling in scramjet engines and found good agreement with experimental results. Thus, the resulting energy equation is implemented in Equation (1).
The governing equations and boundary conditions are shown in Equation (1).
The filtration velocity vector is represented by the three orthogonal velocity components, The asterisk here refers to the dimensional quantities. The terms ν, K, c F , ρ, |u * |, and 'g' in the Darcy-Forchheimer equation are the kinematic fluid viscosity, permeability of the porous material, dimensionless Darcy-Forchheimer form drag coefficient, fluid density, velocity magnitude, and gravity, respectively. The term e z is the unit vector in the z direction. T * , c p , k, β, and λ are the respective temperature, specific heat, thermal conductivity, coefficient of thermal expansion, and a velocity scaling constant. The subscript '0' and 'w' refer to the constant quantity at the top and bottom boundary, respectively. Here, all the properties are assumed constant except the density associated with buoyant forces.
In many transpiration cooling applications, the fuel is often repurposed as the coolant. These fuels, e.g., hydrogen, methane, and other heavier alkanes, are delivered as liquids at high pressures at or above the thermodynamic critical point of the fluid. As the fluid is transpired through the porous material, heating occurs which often causes the fluid to become supercritical (if it is not supercritical to begin with) and causes a reduction in the fluid density. Density gradients are highly nonlinear near the thermodynamic critical point and in the supercritical regions. Accordingly, the stability behavior can be significantly different from systems with linear dependencies. In order to capture these nonlinear effects, density is modeled as a third order polynomial, which qualitatively captures the real-fluid behavior of the coolant as illustrated in Figure 2. This is often considered as the nonlinear Boussinesq approximation [22,23]. Figure 2 shows non-dimensional reduced density, ρ r = ρ ρ c , versus non-dimensional reduced temperature, T r = T * T c , of methane and hydrogen generated using real-fluid equations of state implemented through REFPROP [24]. The subscript 'c' represents the property at the critical point. Figure 2 shows methane and hydrogen for sweep of reduced pressures from P r = 1 to P r =5 -a typical pressure range for coolant delivery [25]. A least squares fit using a third order polynomial is shown in Figure 2 as the solid black line. Figure 3 shows typical conditions cooling the throat of a rocket engine using methane with an at in inlet temperature of 105 K and a wall temperature of 755 K with the density (red-dashed line), the polynomial curve fit (black solid line), and the linear fit (blue small-dashed line) [25]. The values are also used for the fluid properties in Table 1. Equation (2) shows the third order polynomial: where the non-linear portion can be separated as the function f (T) = c 1 T 2 + c 2 T 3 . In Equation (2), T is the non-dimensional scaling of temperature where T = T * −T 0 ∆T , β is an average coefficient of thermal expansion, and ∆T = T w − T 0 . It is worth noting that β is different for the linear case compared to the non-linear case. As a result, it varies from the linear case, where β = 0.00145 to the non-linear curve fit where β 2 = 0.0057. To appropriately compare to linear and non-linear cases, a scale factor, β/β 2 = 3.97, must be accounted for in the results and will be discussed in the later sections.   The system of equations in Equation (1) are non-dimensionalized by introducing the following scaling variables: The system of equations in Equation (1) are rewritten as: is the Grashof number.

Base State
The basic steady state energy equation subjected to the aforementioned boundary conditions with a constant through-flow given by w b = −λ is given by Equation (5), where C = −λ σγ and the subscript 'b' denotes the base state. The basic solution promotes a thermal boundary layer of the exponential form given in Equation (6): This solution is very similar to the one dimensional solution used to model transpirational cooling experimental data given in Ref. [6], which used constant temperatures for both the top and bottom boundary to match experimental conditions. It is worth noting that the non-linearity in the temperature boundary layer is driven by the velocity magnitude, λ. In the limit of the velocity magnitude approaching zero, the temperature boundary layer becomes linear and approaches the same basic solution described in Ref. [17], where a porous material with a horizontal rather than vertical through-flow undergoes heating at the bottom boundary with a fixed top boundary temperature.
The basic density profile can be solved for by inserting Equation (6) into Equation (2) (shown later in Section 5). Accordingly, non-linearity in the density is caused from two independent sources: (1) The through-flow velocity and its effect on the temperature profile and (2) inherent non-linearity in the equation of state. Together, these effects manifest large density gradients at the top of the domain and small gradients at the bottom. The stability of this configuration is analyzed in the next section.

Linear Stability Analysis
The disturbance of Equation (4) is analyzed via linear stability analysis by redefining the velocity and temperature with a base state superimposed with a perturbation, namely: where is an arbitrary small parameter. Next, Equation (7) is substituted into Equation (4) and linearized. The system of equations is further simplified by utilizing the incompressible relation in conjunction with the vertical component of the double curl of the Darcy-Forchheimer equation. The following system reduces to Equation (8): The product of the function g(T b ) and T is a result of linearizing f where g(T b ) = 2c 1 T b + 3c 2 T 2 b . Next, the method of normal modes is introduced where the disturbances take the form < w , T >=< W(z), Θ(z) > e i(α x x+α y y)+Ωt . The system of equations is finally transformed into the following ODE: The ODE can be further simplified if λξ is small. Considering the transpirational flow rates, permeability, and inlet conditions in Refs. [6,25], λ is around 0.01-0.05. While typically ξ << 1. The values for calculating λ and ξ are given in Table 1. Note that eliminating the terms associated with λξ resembles a typical Darcy law, albeit with a nonlinear density. However, the non-dimensional scaling includes viscosity effects which give a Grashof-type number instead of a Rayleigh-type number. Further assuming the principle of exchange of stabilities is valid, the system reduces to Equation (10). For convenience, let

Power Series
The system in Equation (10) will be solved in three steps. First, by considering the case of λ = 0 and g = 0. Second, λ = 0 and g = 0. Lastly, a non-zero λ and g. For the first two cases where λ = 0, the solution can be obtained using a power series expansion. In order to do this, the system is simplified to one equation in terms of Θ.
Note that the recursive formula for a n+4 is a function of a 3 . In order to determine a 3 , the boundary condition D 2 Θ = 0 at z = 1 is used to solve for a 3 . The truncated power series for D 2 Θ will take on the form: where a 3 can be factored out of φ 2 and explicitly solved. The boundary condition for Θ = 0 at z = 1 is then used to numerically solve for Γ as a function of α. Due to expanding the power series at z = 0, high accuracy requires more terms. Truncating the power series to 30 terms results in a maximum percent error of 1.92 × 10 −13 as compared to the analytical solution for Case 1 for wavenumbers from 1 to 6.

Chebyshev Method
Next, a spectral expansion in Chebyshev polynomials is used to solve the system of ODEs given in Equation (10). A detailed description of the method is provided in Ref. [26], and thus, a brief overview is provide here. Employing the Chebyshev technique requires shifting the domain to the Chebyshev domain (−1, 1) and expressing the variables as a finite expansion of Chebyshev polynomials: A weighted inner product with T k is applied to each equation and converted to a generalized eigenvalue problem.
where D 2 are I are the Chebyshev matrix representation of D 2 and identity matrix, respectively. The final eigenvalue system can be solve by passing the matrices through a standard eigenvalue solver, e.g., Matlab's "eig" subroutine.
Determining the onset of instability requires the smallest eigenvalue. Thus, the wavenumber is swept and a Γ determined for every wavenumber. One complexity is added when λ = 0, that is, Γ cannot be separated explicitly as an eigenvalue. In order to overcome this, a 'test' Γ is supplied to the terms containing Γ in the matrix and a 'new' Γ is solved. This iterative process is done until the 'test' Γ and the 'new' Γ converge within a designated value (absolute error <0.001).

Finite Difference Method
Similarly, finite difference method solves Equation (11) as an eigenvalue problem. Since the case of marginal stability is a steady-state problem, the governing equation is classified as an elliptical differential equation. As such, it must be solved by utilizing an implicit method. A system of differential equations is posed, with the governing eigenvalue equation solved at each grid location. For example, Case 1 yields the following: where [D 2 − α 2 2 is the finite difference operator which uses central difference schemes for the various derivatives (D, D 2 , D 3 , D 4 ) before summing them together (D 4 − 4D 3 α + 6D 2 α 2 − 4Dα 3 + α 4 ) over the mesh stencil for each equation. The algorithm proceeds as follows: (1) Prescribe the wavenumber, α.
(2) Fill out operator matrix [D 2 − α 2 2 including boundary conditions listed above. (3) Utilize Matlab's existing "eig" feature to determine possible eigenvalues for a given operator. (4) Find the lowest positive eigenvalue for Γ. For Case 3, the presence of Γ in the operator requires an iterative guess-and-average procedure similar to that described above for the Chebyshev method.

Results
5.1. Case 1: g = 0 and λ = 0 The first case to consider regarding the system given in Equation (10) is the simplest case of λ = 0 and g = 0. Note as λ − → 0 the base state temperature gradient, dT b dz − → −1 ,i.e., the conductive limit. For large λ, the base state approaches that given by Wooding [19]. It is easily demonstrated that the system reduces to Equation (13): It is worth noting that this system is the system given for the Horton-Rogers-Lapwood problem in [14], albeit a Grashof-type number is used (Γ) instead of a Rayleigh-type number.
The analytical solution is simple: It is clear that Γ has a minimum when j = 1, leading to critical parameters of Γ c = 4π 2 and α c = π. Note that Γ includes the effects of σ, which dictates how the porous medium properties ((ρc p ) s ) differ from the fluid properties.
This result can be compared to the three techniques used in in this work and is shown in Figure 4. The signature marginal stability curve is reproduced with all three methods. Table 2 shows the relative error for the three methods compared to the analytical solution between wavenumbers 1 and 6. Thirty terms are used for both the power series method and the Chebyshev method. The power series is the most accurate but has the longest run time. The finite difference method can be improved by using finer mesh instead of the relatively coarse mesh of 100 points, however, the 100 grid points gave comparable run times to the Chebyshev method and thus was selected as a baseline.
Next, the effects of non-linear density are added and the resulting system is shown in Equation (14). As mentioned previously, when a non-linear base density distribution is used, β is different than in Case 1. Thus, a scaling factor is needed to equate the Γ from case 1 to the Γ in Case 2. The Γ for Case 2 will be referred to as Γ 2 and is related to Case one by Γ = Γ 2 β β 2 . The above system is solved using the three numerical methods that were validated for Case 1 above. As shown in Figure 4, these methods again are all in agreement in their prediction of the marginal stability curve. Due to the additional terms in Equation (14), the power series method converges slower. The Chebyshev method converges the fastest out of the three methods and thus is used for comparison. To demonstrate the agreement between the finite difference and Chebyshev, the mesh is increase by a factor of 10. Table 3 demonstrates the degree of convergence between the three methods for Case 2.  Figure 4 shows the marginal stability curve shift to a higher Γ compared to Case 1, indicating a stabilizing effect. Figure 5 helps provide insight to why the non-linear density profile results in a stabilizing effect. Accounting for non-linear changes in density yields a rapid decrease in the initial density near z = 1. Thus, there is a lower amount of heavy fluid over-top the lighter fluid at the bottom (z = 0). The end results is a higher Grashof number occurring at a larger wavenumber. The Γ c and α c for Case 2 are 47.51 and 3.57, respectively. Next, the normalized velocities for each case are compared in Figure 6. The non-linear density base profile causes a shift in the sinusoidal shape toward the top of the domain. The overturning convective velocity currents appears to be biased to where the steepest density gradient occur (see Figure 5). Thus, most of the convective mixing inside the porous media will be shifted towards the top of the domain. In addition, the greater wavenumbers in Case 2 promote tighter convective cell formation. These tighter packed cells will ultimately enhance mixing near the coolant side of the porous media once instability ensues. Lastly, Case 3 is investigated given by the system in Equation (10). The solution method for Case 3 provided here is the Chebyshev method due to the robustness and speed of convergence. It is worth noting that two additional complexities are added when λ = 0. (1) The base state temperature becomes an exponential function and (2) the base state temperature gradient, dT b dz , is also an exponential function. As mentioned previously, the iterative approach described in the Solution Techniques section is implemented. Figure 7 shows the marginal stability curves at various λ's. Typical transpirational flow rates result in a λ ≈ 0.05 [6,25]. Thus, λ is varied from 0.005 to 0.05. Interestingly, increasing λ initially yields a destabilizing effect followed by a stabilizing effect. This is consistent with what was observed Sutton [20] who solved a similar problem to the one in this work (variances include a finite horizontal domain, a linear density Boussinesq approximation, and a different non-dimensionalization). Sutton determined at λ = 0.005 that the onset of instability occurred at 39.42 for a domain with a width twice that of the depth. Here, the minimum point for marginal stability for λ = 0.005 occurs at Γ 2 = 196.49 which corresponds to a Γ = 43.24. Again, as compared to 39.42, the non-linear density provides a stabilizing effect. Figure 7 also illustrates a shift towards higher wavenumbers, which occurs for increasing λ. The critical values for the different flow rates are summarized in Table 4. Table 4 clearly shows the effect of increasing the net through-flow, which causes a consistent shift and increase in the wavenumber resulting in tighter convective cells.
Lastly, all cases are compared in Figure 8. The effect of each constituent complexity can be compared to the original HRL solution. First, the non-linear density dependence increases the stability of the system from a critical value of 4π 2 to 47.51. Then the system is destabilized from 47.51 to 45.73 by adding a base through-flow of λ = 0.05. The system is however regaining stability by increasing the through-flow velocity as shown in Table 4. An increase in wavenumber is observed from Case 1, 2, to 3.

Conclusions
Linear stability analysis was performed for a simplified system modeling transpirational cooling applications. Analysis was performed utilizing three models of increasing detail: Case 1, a Horton-Rogers-Lapwood (HRL) type configuration, Case 2, a HRL configuration with non-linear (cubic) density dependence on temperature, and Case 3, both a cubic density function and a base through-flow velocity. In all cases, linear stability analysis was performed to yield governing equations for the system. These equations were then solved utilizing multiple numerical methods, comparing results from a combination of analytical solutions, finite difference, power series, and Chebyshev methods. Results demonstrated excellent consistency in predictions across these methods.
The addition of a cubic density dependence on temperature had several effects as compared to Case 1 (HRL problem). First, accounting for the non-linear density variation during the transpirational heating process shifts the marginal stability curve to a higher Grashof-type number resulting in increased stability. This is due to sharp density gradients near the top wall decreasing the density such that an overlaying fluid is not heavy as compared to the linear density case. In addition, a shift to higher wavenumbers is observed. The cubic density profile also shifted velocity peak toward the top of the porous media. Finally, Case 3 studied the effects of the addition of a base through-flow velocity-the base flow is in the same direction as the gravity vector. A base through-flow velocity sweep revealed that initially increasing the base through-flow destabilizes followed by a re-stabilizing effect. However, it is worth noting that all base flows for Case 3 considered in this work were more unstable as compared to Case 2, i.e., the highest base velocity resulted in a lower critical Grashof-type number as compared to Case 2. In both Cases 2 and 3, accounting for the non-linear variation in density increases the stability of the system.
The implications of an additional mode of convection (buoyancy-induced convection) is increased mixing, i.e., increased heat transfer inside the porous medium. Accounting for this additional mode of convection in cases where instability ensues is important when it comes to thermal management of the system. Thus, an awareness that sharp density gradients at the top of the domain from non-linear density variation in the system, as well as increased through-flow have a stabilizing effect on the system are pertinent to design and modeling considerations.
Author Contributions: The framing of the problem was carried out by C.T.W. The governing equations, non-dimensionalization, perturbed equations, and base state solutions were independently generated by both authors. C.T.W. implemented the power series and Chebyshev solution techniques while P.R.J. developed the finite difference scheme and solution. The thermodynamic property analysis using real-fluid equations of state along with the non-linear density implementation was completed by C.T.W. Both authors contributed to the writing, editing, and figure generation of the manuscript. All authors have read and agreed to the published version of the manuscript.