Abstract
This paper considers the numerical optimization of a double ramp scramjet inlet using magnetohydrodynamic (MHD) effects together with inlet ramp angle changes. The parameter being optimized is the mass capture at the throat of the inlet, such that spillage effects for less than design Mach numbers are reduced. The control parameters for the optimization include the MHD effects in conjunction with ramp angle changes. To enhance the MHD effects different ionization scenarios depending upon the alignment of the magnetic field are considered. The flow solution is based on the Advection Upstream Splitting Method (AUSM) that accounts for the MHD source terms as well. A numerical Broyden-Flecher-Goldfarb-Shanno- (BFGS-) based procedure is utilized to optimize the inlet mass capture. Numerical validation results compared to published results in the literature as well as the outcome of the optimization procedure are summarized to illustrate the efficacy of the approach.
1. Introduction
Scramjet engine inlet flow is subject to many engineering tribulations, mainly due to the fact that the base geometry of the engine is suited to a very narrow range of flight conditions. In order to have the engine operate efficiently across a broader range of flight conditions the inlet flow must be adjusted. One obvious method to tune the inlet flow is a mechanically actuated, variable geometry inlet that can adjust its shape in flight to achieve optimal inlet conditions [1–3]. Since there are pressure losses across the inlet, the shape of the inlet needs to be optimized to minimize these losses. Another novel method of inlet flow optimization involves the use of magnetohydrodynamics (MHD) to control the incoming flow instead of the mechanically actuated surfaces.
The scramjet engine is designed for hypersonic flight with supersonic combustion and is flown at speeds ranging approximately from Mach to . At these speeds the air ahead of the vehicle can become ionized thereby making a case for magnetohydrodynamic flow control. Since the ionization levels in this type of flow are quite low, it is considered in the MHD community as a “cold flow.’’ Several methods have been explored in increasing the ionization of this type of “cold” flow including seeding and electron beams. Work done by Macheret and Miles suggests that the most efficient method is the electron beam [4–6]. In the limits of increasing conductivity in the flow, increasing the magnetic field strength is of course also very restricted.
It is to be mentioned that initial studies of optimizing duct flows for MHD power generators [7, 8] later evolved into the “AJAX’’ concept proposed by Russian scientists as a means to provide heat protection for hypersonic flight vehicles [9]. The resulting enthalpy reduction due to an MHD generator in the inlet of the engine could theoretically allow the operation of a conventional turbojet or ramjet engine in hypersonic flight instead of a scramjet engine. If the enthalpy extraction was not sufficient to avoid the need for a scramjet engine, the MHD power generator could still lower inlet temperatures to tolerable levels. As well the addition of a bypass system to further accelerate the exhaust gases could increase the specific thrust of the engine during off-design conditions and flow control [10–14].
Previous studies on the effects of MHD control on inlet flow optimization [15] were conducted on a symmetrical double ramp inlet as depicted in Figure 1. These studies involved two cases, first with the magnetic field oriented -plane, or into the page, and secondly with the magnetic field implemented in the - plane. A number of simulations were run with varying conductivity and magnetic field angle and in all cases no improvement was found in the pressure recovery. The application of a magnetic field to the charged flow always results in an increase in the pressure loss. This is consistent with the theory in that energy is not being added to the flow and the resultant Joule heating is only a detriment to the pressure recovery.
The main purpose of this study is to increase the mass capture of a double ramp cowl style scramjet inlet. The performance metric is optimized using a combination of changing geometry (ramp angles) and including MHD effects. The discretized Euler's equations are augmented with MHD source terms and the flow solution is obtained via the AUSM method. The flow solution methodology is validated against test cases drawn from literature available in the public domain. The optimization methodology is also validated together with the flow solver on the traditional double ramp inlet problem with and without MHD terms. Finally, the optimization results for the combined MHD inlet geometry are presented.
2. The Cowl Style Scramjet Inlet
The double ramp scramjet inlet studied in this paper is a cowl style inlet. Such an inlet represents a forebody external compression region followed by a small region of internal compression. This mixed compression configuration can balance the problems of high external drag in the case of full external compression and excessive viscous effects during full internal compression. The optimal cowl inlet configuration is the well-known “shock-on-lip’’ condition shown in Figure 2.
During off-nominal flight conditions when the flow Mach number is less than the design value, the shocks will move ahead of the cowl lip and some of the compressed air will escape the inlet resulting in “spillage’’ and a decrease in the mass capture. In flight conditions where the flow Mach number is greater than the design value, the shocks move into the inlet causing multiple reflected shocks, loss of total pressure, possible boundary layer separation, and engine unstart [16]. It is proposed that the optimal “shock-on-lip’’ configuration can be recovered via flow control employing magnetohydrodynamic source terms as in [16]. As such, the flow control methodology in this paper compares two styles of ionization. The first is inspired by the work from Shneider et al. [16] who focused on a “virtual cowl’’ scenario in which localized off-body energy addition was used to increase the mass capture and pressure recovery at Mach numbers less than the design value. In their work the virtual cowl is a heated region upstream of and slightly below the cowl lip. This heated region acts as a scoop by deflecting the incoming flow and increasing mass capture [17].
The second ionization scenario was presented by Shneider et al. [16, 18] as well as Sheikin and Kuranov [19] who implemented a uniform magnetic field in which they assume an ionized region enclosed by lines that are parallel to the magnetic field lines to improve scramjet inlet performance at off-design conditions. Figure 3 demonstrates these two scenarios. The first scenario is characterized by the stationary-ionized region upstream of and centered on the cowl lip. The magnetic field is applied to this charged region with various angles to determine the capabilities of influencing the inlet mass capture. The second scenario is characterized by a moving ionization region which is coincident with the magnetic field implementation. In this case the magnetic field and ionized regions will be applied together at various angles in an attempt to influence the flow.
(a) Moving magnetic field and stationary ionized zone
(b) Moving magnetic field and coincident ionizing e-beam
Also, it has been independently shown in works such as in [1, 2] that geometric optimization of inlet shapes (namely, the ramp angles in this case) could lead to optimal pressure recovery as well as mass capture. This paper investigates the possibility of simultaneously optimizing the geometry of the inlet (as was done in [1, 2] and the ionization beam angle [16, 18, 19]). The geometric optimization portion of the work will adjust the two ramp angles, and .
3. Euler Equations with MHD Effects
The equations of motion that govern inviscid, compressible fluid flow in a region are given by Euler's equations [20] and are summarized as follows:
These equations can be recast into the following form in order to facilitate the implementation of a flux splitting flow solver method used in this study. The equations are limited to 2D and are consistent with the formulations employed for the study of inlets: where is called the flow solution vector, and are known as the flux vectors, and represents the source terms: In addition, the following perfect gas relations are assumed to hold
In modeling the source terms we consider the effects of a charged flow through the inlet with an applied electromagnetic field. This can be accomplished with the addition of appropriate electrodynamic terms [21] to Euler's equations. It is to be mentioned that the MHD flows through the inlet are characterized as having a very low (0.0001 1) magnetic Reynolds number, . (. High magnetic Reynolds numbers (1) are characteristic of fusion research and astrophysical phenomena.)
This suggests that the conductivity in MHD flows is very low and therefore the current and hence the induced electric field are also very small. This allows us to assume to be constant; therefore
As a consequence of the above we may introduce a scalar electric potential such that and . The current density is then calculated using Ohm's Law: Thus we can then calculate the electric potential for a given magnetic field and flow conductivity as follows [22–24]: Therefore, where is the elemental area. This allows us to then calculate the current density directly from (6), which will become necessary as we implement the source terms into the equations. Note that we neglect the effects of ion slip in this model.
With the application of a magnetic field to a charged flow the body forces and volumetric heating effects are no longer negligible. The body force term known as the Lorentz force is given by the vector , and the volumetric heating known as Joule heating is given as . As seen earlier in the development of the Euler equations the contribution to the momentum equation is strictly due to the Lorentz forces where the contribution to the energy equation is the total rate of energy addition, due to both volumetric heating and work done by the Lorentz forces.
In (2) the source term, is now given as
The left-hand side of (2) is unchanged and is hyperbolic for Mach numbers greater than one. The right-hand side however is unconditionally elliptic for smooth variations of material properties; we therefore need to implement a Poisson solver to obtain the source terms. If however we implement a magnetic field in the - plane and consider an ideally sectioned Faraday MHD generator such that the Hall effect is neutralized, the resultant current density is only in the -direction. This greatly simplifies the computations as well as the implementation and is used in this study to investigate the feasibility of simultaneous MHD and geometric optimization: Therefore, For this study we set to zero which corresponds to a short circuit of the electric field.
Thus, the source terms are then trivially found as follows without the need for a Poisson solver:
4. Flow Solution
4.1. Numerical Approach
We employ the flux splitting method developed by Liou and Steffen [25], known as the advection upstream splitting method (AUSM). It is a first order scheme that is relatively simple to implement and yet still has the ability to resolve shock structures. It requires only operations in contrast to operations needed for a Roe splitting scheme (where is the number of equations). The choice of the scheme is motivated by the simplicity and ease of implementation for a feasibility type study. Of course more accurate solutions can be obtained using higher-order schemes. The scheme is developed as follows: the flux vector is split into two components, convective and pressure terms: The convective terms are propagated at the cell interfaces by an appropriately defined velocity and the pressure term is propagated at acoustic wave speeds. This leads to the two terms being discretized separately. For the convective terms at an interface , where The Mach number splitting method for the left and right states utilizes Van Leers definitions as follows:
Similarly for the pressure terms, where the pressure splitting is weighted using the second-order polynomial of the characteristic speeds as This splitting of the advection and pressure terms allows for the complete definition of the inviscid flux vector.
4.2. Grid Generation
The grid chosen for this study is a simple algebraic style grid with a Thomas Middlecoff control function applied for smoothing. The algebraic grid utilizes uniformly spaced grid points in the -direction and a transformation in the -direction which allows for clustering of the grid points near the wall boundary [26]: As the stretching parameter approaches , this transformation clusters more points near . In order to implement this transformation we need the inverse of the above equations given as
4.3. Boundary Conditions
In setting the boundary conditions for these simulations we must consider that the equations are hyperbolic in nature for Mach numbers greater than one. We therefore set the inlet conditions equal to those of the free stream and at the exit we use conditions from the interior. Wall boundary conditions on the upper and lower surfaces are reflecting surfaces where the flow normal to the surface is mirrored.
5. Optimization
As mentioned earlier, in this study the design variables are the ionization beam angle that characterizes the MHD effects and the ramp angles , that represent the geometry of the inlet. The objective function that is being optimized is the mass capture given by where the , indices refer to the mesh points on the exit plane of the inlet. The constraints imposed in this optimization procedure are the flow equations, solved by use of the flow solver described earlier.
Due to the high computing time of the optimization code we first ran a few broad parametric studies to establish the trends of each ionization method followed by the geometry changes. Once we had a general region for an optimal configuration we implemented the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method in the form of the “fmincon’’ function in MATLAB (MATLAB/SIMULINK is the trademark of The MathWorks Inc; Natick, MA,USA) with initial conditions close to the values obtained via the parametric studies. Note this step could be automated via a genetic algorithm-like approach that characterizes the feasible domain of parameters. The function “fmincon’’ is designed to find the minimum of a constrained nonlinear multivariable function [27]. It uses a Sequential Quadratic Programming (SQP) method where an estimate of the Hessian of the Lagrangian is updated at each iteration using the BFGS formula. The Lagrangian is the cost function augmented with the constraint function as shown below:
In our study, the function “fmincon’’ was found to be very sensitive to the step size as well as the scaling of the objective function itself. Several combinations of the above mentioned variables were tried in order to obtain the best results and it was found that a scaling of provided good convergence properties.
The implementation of this optimization routine is laid out in Figure 14. As is shown in the flow chart the flow solution is calculated times for each update of the Hessian, where is the number of design parameters. The values of the objective function are stored from iteration to iteration for determination of the minimum and calculation of the gradients.
5.1. Validation of the AUSM Scheme: BFGS Scheme-Based Optimal Pressure Recovery for the Double Ramp Shock Canceled Inlet
Validation studies included validating the AUSM scheme and comparing the geometric optimization of the double ramp inlet to results published in the literature. A symmetrical two ramp inlet, represented in Figure 4, was used in previous studies by several research groups as a base model for optimization procedures. Korte and Auslender [1] used it to demonstrate a new optimization procedure utilizing a least-squares parabolized Navier-Stokes solver. This style inlet was also utilized by Munipalli et al. [2] for validation of their CFD and optimizer code. Both groups then went on to optimize a contoured inlet defined with cubic splines.
(a) Initial inlet configuration with shock canceling
(b) Mach contours for the shock canceled configuration
In previous work by the authors [15], the inlet was used to demonstrate the use of numerical gradients in optimization of the geometry of the inlet during off-design conditions as well as the preliminary investigations into the possible uses of MHD control for flow optimization.
This symmetric inlet was optimized for the mass-weighted total pressure loss from free stream to the throat. The optimal geometry of the inlet allows for the shocks to reflect onto the opposing side of the throat. For operation at speed greater or less than design speed the pressure losses will increase. The initial configuration is the shock-canceled case where the shocks reflect at the centerline and converge on to the throat. The initial ramp angles do not represent an optimal configuration for any Mach number; they are simply a nonoptimal case for our given Mach number.
Figure 4 demonstrates the initial conditions for the inlet and Table 1 gives the analytical solution for the initial state of the shock canceled inlet. From this table we can see that total pressure recovery of the shock-canceled inlet is and will be used in validation of the CFD flow solver.
The early version of the code was validated in our previous work [15] and redone here using the total pressure recovery of the shock-canceled condition for the two-ramp inlet. This case as demonstrated in Table 1 is determined analytically based on the pressure recovery factors going from one region to the other across the inlet. The corresponding angles degrees and degrees and a total pressure recovery of are thereby obtained. The resulting flow solution described here gives a maximum total pressure recovery of for the above angles, as compared to the result of from Munipalli [2] and from Korte [1]. This is a reasonable result when considering the first order accuracy of the flow solver and the much simplified grid used in this case.
We note that the results of this validation, given in Table 2 are consistent with those of Korte and Auslender. The initial and final configurations are summarized in Table 2. Figure 5 is a plot of the Mach number contours for the final optimized configuration.
6. Optimization Results for Simultaneous Ionization Beam Angle and Geometry Changes Applied to the Cowl Style Inlet
Buoyed by the confidence in the optimization procedure and the flow solution after validation against published results (see earlier section), the procedure was applied to simultaneous optimization of the ionization beam angle and the ramp geometry. This section summarizes the results obtained for the Cowl style inlet (the main focus of the study). The results of these trials are shown in Figure 6 where deg, deg, and the mass capture equal to kg/s.
To simulate a less than design Mach number flow we adjust the ramp angles to deg and deg which results in a mass capture of kg/s. See Figure 7 for the Mach contours in this case.
Given this off-nominal design condition, we investigate the ability of an applied magnetic field to direct the flow back to the optimal mass capture configuration. Two different scenarios are considered: () a moving -beam type ionization method [19] such that the magnetic field and ionization region are movable and coincident and () a moving magnetic field but stationary-ionized region [16] (Figure 3 demonstrates each scenario).
Figures 8, 9, 10, and 11 show the ability of an applied magnetic field to direct the flow back to the optimal mass capture configuration. It is clear from these results that the larger the magnetic field strength and the larger the conductivity, the greater the influence on the flow. It is interesting to note that the stationary-ionized zone has a much greater ability to manipulate the flow than the coincident-ionized region. For the largest values of and (line corresponding to in Figure 11), we can see that the mass capture is indeed approaching that of the optimal situation.
(a) Coincident e-beam ionization
(b) Stationary ionization
(a) Coincident e-beam ionization
(b) Stationary ionization
(a) Coincident e-beam ionization
(b) Stationary ionization
(a) Coincident e-beam ionization
(b) Stationary ionization
As mentioned before, based on the results from a broad parametric study, the problem and the optimization procedure was set up. To account for the sensitivity of the optimization routine several different initial conditions were evaluated to see how well the results converged. Finally, the initial conditions for the angle were chosen above and below the approximate optimal values given by the parametric study. Tables 3 and 4 summarize the results obtained from these numerical experiments.
Comparing the results of the optimizer with those of the parametric study for the stationary ionized zone we can see a very nice correlation. The spread in the optimized results is not large and is consistent with the parametric study. However in comparing the optimizer results with those of the parametric study in the case of the coincident ionization zone we see a bit of discrepancy. There is a large spread in the values given by the optimizer depending on whether we began the optimization above or below the peak value shown in the parametric study.
Additional Remarks
It is to be mentioned that a limitation of the current study is that it ignores the specific effects of power deposition due to the electron beam. The optimization problem was solved with the intent of finding a solution for the most optimistic scenario and as a result the power deposition was ignored to keep the simulation complexity low. To model the aspect of power deposition [28], it would require one to include a more sophisticated nonequilibrium temperature model that includes distribution of energy due to particle collisions between the e-beam particles and the gas molecules. Additionally, delays in the energy distribution need to be accounted for. It is also to be mentioned that the power deposition due to the e-beam is a function of the location where the e-beam is applied and the local pressure conditions. As such, while the incorporation of the power deposition definitely strengthens the overall analysis, it does not take away from the fact that there is some optimization achieved within a very optimistic scenario of reduced energy losses vis-a-vis the neglected energy losses due to the power deposition.
In order to better understand this result we conducted another parametric study around the peak value. We limited the study to the Tesla and mho/m case and varied the magnetic field angle from to degrees. As can be seen in Figure 12 the results of this study show that the curve is not very smooth and this in some measure explains the large range in values given by the optimizer for this case.
Finally, Figure 13 shows the outcome of a simultaneous optimization of the geometry and the ionization beam angle. In reality, while this would require actuating the ramp surfaces as well as an ionization method and a way to generate the magnetic field, it is also represents the possibility of fine tuning the magnetic field. Again the off-nominal case of deg and deg was used with the initial magnetic field angle of deg. We limited this study to the stationary ionization zone with a conductivity of mho/m and a magnetic field strength of Tesla.
(a) Magnetic field angle history
(b) Geometry angle history
Figure 13 shows the design history of the optimization routine at each iteration of the inside loop in Figure 14. The outside loop iterations or number of times that the Hessian was updated is equal to the total number of flow solver iterations divided by the number of design variables plus (). From Figure 13 we can tell that there were approximately iterations of the flow solver and design parameters, giving us approximately updates of the Hessian. The values plotted are scaled as per the scaling function mentioned earlier.
The final configuration for the geometry is very close to that of the optimal mass capture with no MHD source term present (see Figure 6), and the final magnetic field angle is consistent with that of the previous case providing a high level of confidence in our procedure and implementation. The results are also summarized in Table 5.
7. Summary and Conclusions
This paper studied the numerical optimization of a double ramp cowl style scramjet inlet using magnetohydrodynamic (MHD) effects together with inlet ramp angle changes. An AUSM-based flow solver was utilized to solve the 2D inviscid, compressible Euler equations subject to MHD source terms. The objective function in the optimization was the mass capture at the throat of a cowl style inlet, so that spillage effects for less than design Mach numbers are reduced. The optimization procedure implemented in this study was a numerical Broyden-Flecher-Goldfarb-Shanno- (BFGS-) based procedure. Numerical validation results compared to published results in the literature have been summarized at various stages that include flow solution and the optimization procedure. It is shown that spillage occurring from off-nominal geometries can be reduced by employing MHD control. We also demonstrate a more attractive case of spillage reduction employing simultaneous optimization of the ionization beam angle and the ramp angles.
Nomenclature
| & : | Ramp angles (degrees) |
| : | Magnetic field strength (Tesla) |
| : | Direction of the e-beam measured from -axis. (degrees) |
| : | Fluid density () |
| Fluid momenta in , , directions () | |
| : | pressure () |
| : | Energy (Joules) |
| : | Velocity () |
| : | External force (N) |
| : | Flow solution vector |
| : | Flux vectors |
| : | Source term |
| : | Ratio of specific heats |
| : | Temperature () |
| : | Gas constant () |
| : | Speed of sound () |
| : | Magnetic Reynolds number |
| : | Permeability of free space () |
| : | Gas conductivity () |
| : | Characteristic length (m) |
| : | Electric field intensity () |
| : | Electric potential () |
| : | Current density |
| : | Total pressure at the th station |
| : | Objective function |
| : | Lagrange multipliers |
| : | Constraint function. |
Acknowledgment
The authors gratefully acknowledge HyperComp Inc. (POC. Dr. Ramakanth Munipalli) for the financial and technical support for this work.