Abstract
The objective of this paper is to numerically study the shock wave/boundary layer interaction and boundary layer separation. The first stage of this research is the development of methodology, flow simulations, and data analysis. When comparing the plots, it can be seen that the results of the check of the methodology were similar. Following, methodologies were developed and simulations were carried out considering the compression corner model. It was noticed that the shock wave could be identified by the jump on the pressure profile near the leading edge and by analyzing the thermodynamic properties of the plate. An increase in pressure, flow inversion, and boundary layer separation through negative values of the friction coefficient was observed, and negative speed at the wall was observed due to the presence of a plateau on the pressure curves. Flow expansion and further reattachment of the boundary layer were also seen. It is possible to observe type VI shock-shock interference and the triple point T, causing a series of expansion waves to form. Finally, an increase in the Mach number, a decrease in the corner compression angle, and a decrease in the wall temperature interfere and reduce the possibility of separating the boundary layer.
1. Introduction
When a vehicle is at supersonic or hypersonic speeds, there is the formation of shock waves and their interaction with the boundary layer. The shock wave/boundary layer interaction has been studied for the last 80 years, due to the advancement of space programs, its presence in supersonic flights, and the need to solve structural and heat transfer challenges. The first research in the 40s showed that there is a complex flow field during the boundary layer interaction with the shock wave.
Dolling [1] carried out a bibliographical review of the research from the previous 50 years on the shock wave/boundary layer interaction. Regarding heating, the author already reported that the heating related to the shock wave/boundary layer interaction can be 10–100 times greater than under the incoming attached boundary layer flow. It has been observed that this subject has made a lot of computational advancement during these years but still has challenges of accuracy in predictions. The author hoped that the computational advances in the next 10–20 years would resolve many issues of concern. Indeed, the authors were right, and between 2020 and 2023, there were considerable computational advancements and solution methods for flows.
The interaction of a boundary layer with a shock wave produces several flow features and should not be overlooked for aerospace vehicle design. The problems involve external and internal flows in aircraft and rockets, which can cause control loss, peaks of thermomechanical loads, and instabilities.
The separation of the boundary layer can cause unsteadiness in the flow, generating oscillations that can damage the engine. This unsteadiness involves very low frequencies and can occur at high frequencies in turbulent flows and/or in separation bubbles. In totally supersonic cases, instabilities generally do not occur, as the external flow tends to isolate from downstream disturbances [2]. For unsteadiness cases, control methods have been studied.
The maximum mean and fluctuating pressure levels and thermal loads are generally found in regions of shock/boundary layer and shock/shear layer interactions and can affect vehicle and component geometry, structural integrity, material selection, fatigue life, the design of thermal protection systems, weight, and cost [1]. In internal flows, as in the scramjet engine, it can increase pressure losses and distortion, inversion, and sudden changes in the flow that can cause unstart of the air inlet [3]. The boundary layer can decrease the inlet mass flow rate which reduces the amount of oxygen that would be used for combustion, leading to, for instance, loss of thrust, specific impulse, and even vehicle control.
Figures 1 and 2 represent the behavior of supersonic flow over a ramp without and with boundary layer separation, respectively. In Figure 1, it can be seen that the pressure increase due to the interaction of the boundary layer and the ramp shock causes the subsonic layer to expand, forming compression waves () and coalescing into the shock (). In Figure 2, the deflection (α) of the ramp exceeds the force that the boundary layer can withstand, separation occurs at point S, and the formation of the separation shock () occurs. After the point of separation, there is a flow recirculation region, and then, there is the reflow of the flow in R and the re-replacement shock (). Then, the shocks intersect at point T.
Boundary layer separation and its intensity depend on several parameters, for instance, the ramp angle, Mach number, Reynolds number, wall temperature, and pressure gradient. In fact, the prediction of laminar boundary layer separation has been the subject of a great deal of papers. Lange [4] showed experimental results of the Mach number and Reynolds in laminar and turbulent boundary layer separation and compared it with available methods to predict separation. Lewis et al. [5] carried out an investigation of the boundary layer separation boundaries in the compression corner for Mach 6 considering the constant temperature and adiabatic wall. A good approximation for the adiabatic wall was found, but it was not possible to correlate with cooling. Hodge [6] compared the results obtained by the MacCormack method with experimental results. The code obtained good results for cases with no separation and cases with little separation. However, the code had anomalies and did not achieve similar performance for separate and recirculating regions.
This topic has been studied for a long time. More recently, in 2014, Gaitonde [7] carried out a review on the topic and showed that current objectives include reduction of the separation region, surface charges, and modulation of spectral content and that SBLI studies have benefited greatly from scramjet test flights.
The objective of the present work was to numerically investigate the shock wave/boundary layer interaction in the flow over a compression ramp. These conditions were simulated on supersonic Mach numbers and deflection angles considering constant wall temperature and adiabatic wall conditions. This investigation seeks to characterize the physical phenomena that occur during shock wave/boundary layer interactions and understand in which conditions the flow deflection angle, Mach number, and wall temperature could cause boundary layer separation (and separation length).
The text is organized as follows: Section 2 reports the methodology used, where the MacCormack method is applied. Section 3 shows the results and analysis for Mach number, temperature, pressure, skin friction coefficient, heat flux, and numerical schlieren. Finally, Section 4 shows the conclusions.
2. Methodology
2.1. Navier–Stokes Equations
This work used the explicit finite difference solution of the two-dimensional Navier–Stokes equations: continuity, x momentum, y momentum, and energy, neglecting body forces and volumetric heating [8]:where the column vectors U, E, and F were obtained by
The equations are functions of the pressure (), velocity components (), density (), shear stresses (), total energy per unit volume (), and heat fluxes ().
is physically composed of two parcels, one accounts for the internal energy and the other one for the kinetic energy:where is the internal energy per unit mass and .
The components of the heat flux vector and (following Fourier’s law) were calculated as functions of the temperature () and the thermal conductivity ():
Considering relations for Newtonian fluids of Stokes [9], a laminar flow and shear stresses [9] were represented in terms of velocity gradients as follows:
and are the dynamic viscosity and volumetric viscosity coefficient.
The volumetric viscosity coefficient was calculated using
The system of equations consists of four equations with nine unknowns. Therefore, the following considerations were necessary to solve the system: (i) an ideal gas (equation-of-state), (ii) a calorically perfect gas, (iii) Sutherland’s law to calculate dynamic viscosity, and (iv) a constant Prandtl number which is equal to 0.71.
2.2. MacCormack Method
The MacCormack method is a finite difference scheme, in which there is a time walk to get a steady state solution. Thus, the properties are obtained at each spatial location (i, j) of the discretized domain. According to Anderson [8], the MacCormack method can be applied using two steps, a prediction step of time and a correction step with second-order precision at time .
In the prediction step, where forward differences in x and y were applied, the solution vector was obtained by
The solution vector for the correction step was calculated with backward differences in x and y as follows:
The MacCormack method has an explicit formulation, so it was necessary to use a stability criterion for the time step. In this work, a version of the Courant–Friedrichs–Lewy (CFL) condition described by Anderson [8] was used, where the time step () was calculated by
The Courant number works as a correction factor added to the calculations to adapt them to simulation observations. The Courant number’s recommendations vary as the types of simulations vary, but the most recommended practice is using values less than 1. In fact, this work used .
Vector U was decoded after the prediction and correction steps to obtain the primitive variables: density, velocity components, and internal energy. The convergence test was performed comparing the values of density at each point of the domain (). Numerical validation was obtained when the deviation was less than 1% between the input and output mass flow of the computational domain. The mass flow was found by mass conservation.
The computational domain has a rectangular-structured configuration and a length of 0.00001 m, a small value, but sufficient for the desired physical phenomena and large at the molecular level [8]. The program was developed in the MATLAB language [10].
The domain size must be sufficient for the shock wave to be inside, and for the trailing edge, Blasius prediction was used [8]. Boundary layer thickness and vertical height were determined by
2.3. Numerical Analysis
The skin friction coefficient and the Stanton number were calculated by
The Stanton number and the friction coefficient were compared with the methodology of Simeonides [11]. The hypersonic viscous interaction parameter, the reference temperature, and the recovery temperature were defined in equations (14)–(16). The recovery factor (r) for laminar flow was calculated using .
The numerical schlieren [12] was analyzed through the contrast:
2.4. Compression Corner
2.4.1. Coordinate Transformations
The MacCormack method in the previous section was developed on a rectangular domain. Thus, it was necessary to implement the methodology variables into a transformed coordinate system, which is necessary for the study of the flow over the compression corner. In the coordinate transformation, the equations of the correction and prediction steps were adapted for the new coordinates, η and ξ. Likewise, normal and shear stresses and heat flux components were obtained for the new coordinate system.
According to Anderson [8], and can be defined bywhere and were defined by
The solution vector for the prediction step is
The solution vector for the correction step is
The stresses can be written as [13]
Heat flux vectors were written as follows:
Furthermore, using and , schlieren expression can be written in relation to and as follows:
2.4.2. Physical Domain
Figure 3 represents the physical domain for the coordinate transformation according to Anderson’s test case [8]. It was adopted so that the expansion corner can start at the midpoint of the domain.
Case 1 represents the leading edge, where freestream pressure and temperature were used under no-slip conditions . Case 2 represents the top and left boundary conditions (except for the leading edge), and the properties of freestream were used. Case 3 represents the interior points of the exit plane, and the properties were calculated using extrapolation. Case 4 represents the boundary condition of the plate surface (except for the leading edge). The nonslip condition was considered (), and the wall temperature () and pressure were calculated using extrapolation. In the case of the adiabatic wall, .
and can be calculated by equations (19) and (20), while varies from 0 to L and varies from 0 to 1.
Based on Figure 3, one can conclude that when ,so that
On the other hand, when ,
Therefore,
3. Results
3.1. Preliminary Analysis
The first step was to validate the results obtained by the test case described by Anderson [8]. Thus, the simulation of a supersonic flow on a flat plate without the angle of attack and without incident shock waves was performed, using the MacCormack method for the constant temperature and adiabatic wall. The computational domain and freestream properties are described in Tables 1 and 2, which are domain size, angle of attack, CFL number, plate length (), domain vertical height (), step in directions x () and y (), Mach number (), temperature (), pressure (), ratio of specific heats (), Reynolds number (), and Prandtl number ().
Figures 4–7 show the nondimensionalized pressure and temperature distribution along the surface of the flat plate and the trailing edge. When comparing the plots, it can be seen that the results are similar to those obtained by Anderson [8] for constant temperature conditions and adiabatic walls.
The verification of the flow properties was performed in the isothermal wall condition. Figures 8–11 show the distributions of the temperature, Mach number, pressure, and numerical schlieren along the nondimensionalized computational domain in relation to the plate length and the vertical height of the domain .
The flow enters from the left at Mach 4, and the formation of the boundary layer at the leading edge of the plate can be seen. Boundary layer displacement due to viscous effects causes the formation of an oblique shock wave. It can be seen that, on the wall, the temperature is constant and the velocity is zero. The temperature has the highest values near the plate surface due to viscous dissipation of kinetic energy within the boundary layer.
Temperature, Mach number, velocity, and pressure profiles are represented along the nondimensionalized computational domain in relation to the final length and the vertical height of the domain . It can be seen that the profiles capture the shock wave through the jump in the curve. In Figure 12, for example, it can be seen that the height of the shock wave increases as the trailing edge approaches for and for ).
It is observed that close to the plate, due to the formation of the boundary layer, the velocity is null due to the non-slip condition. Thus, velocity has smaller values, and temperature and pressure have larger values near the surface of the plate. Toward the trailing edge, the velocity gradually increases to the freestream condition. Thus, the pressure and temperature also decrease for the freestream values. Analyzing the velocity at the positions of for 0.01449, 0.5, and 1, it is noticed that the velocity reduces toward the trailing edge (). This is because increasing the thickness of the boundary layer causes a reduction in the velocity gradient (see Figures 12–15).
The profiles of the Stanton number, skin friction coefficient, and heat flux along the nondimensionalized computational domain are represented in Figures 16–18. These properties are greatest at the leading edge and decrease at the trailing edge. This happens due to the passage of the oblique shock wave and the increase in temperature.
3.2. Code Validation
Code validation was performed by comparing the results of the coefficient of friction (Figure 19) and the Stanton number (Figure 20) along the plate and changing the Mach number for (Figures 21 and 22) with the semiempirical methodology of experiments of Simeonides [11] for Mach numbers 4 and 7. The properties of the computational domain input conditions used in this validation are shown in Table 3, with the constant wall temperature.
Near the end plate, the values come closer in both methodologies. It was noticed that the results tend to approach Mach 5. From Mach 5 onwards, as the Mach number increases, the values of the Stanton number and the skin friction coefficient tend to move away from the two methodologies (0 deg). Comparing the simulation data with the experimental data, a maximum difference of 17% of the skin friction coefficient (Figure 21) and that of 6.9% of the Stanton number (Figure 22) were found, and these errors indicate good agreement.
Table 2 shows the number of iterations necessary for convergence for different sizes of the computational domain. It is observed that, in the 100 × 100 domain, there is an increase of 98.19% in the number of iterations in relation to the 50 × 50 case. In the 150 × 100 domain, there is an increase of 46.82% in relation to the previous case, and in the 200 × 200 domain, there is an increase of 30.9%. Thus, the number of iterations tends to converge with an increase in the computational domain Table 4 shows the number of iterations necessary.
The graph in Figure 23 shows the relationship between the domain size and the ratio between the maximum temperature in the computational domain and the freestream temperature (). It is observed that the temperature ratio in the 50 × 50 domain is 97.42% of the 100 × 100 case; in the 100 × 100 domain, it is 97.87% of the 150 × 150 case. In the 150 × 150 domain, it is 99.21% of the 200 × 200 case.
3.3. Compression Corner
The simulation of supersonic flows was applied to the compression corner case, with the test matrix depicted in Table 5. At the constant temperature ) and adiabatic wall, the Mach number varied from 3 to 6 and flow deflection angles varied from 15 to 30 deg.
Figures 24–37 show the temperature, pressure, wall pressure, Mach number, skin friction coefficient, Stanton number, and heat flux with the constant wall temperature and adiabatic wall. The x-axis coordinate was nondimensionalized with respect to the total flat plate length (), and the y-axis coordinate was nondimensionalized with respect to the vertical height of the domain .
The temperature (Figures 24 and 25) decreases along the boundary layer, and it has a small growth curve as it passes through the interference shock wave and then converges to the freestream value. The pressure (Figures 26 and 27) remains constant at the end of the plate within the boundary layer, has abrupt decay on passing through the shock, and converts to freestream pressure. The same effect occurs with the Mach number (Figures 30 and 31), where the fluid accelerates as it approaches the end of the boundary layer, and it has decay when passing through the shock and then reaches the freestream Mach number value. The cases with the adiabatic wall condition reach maximum values of temperature and pressure much higher than the cases with constant temperature, about 2000 K and 1.8 MPa against 700 K and 0.8 MPa at the end of the ramp.
Near the plate’s leading edge, the wall pressure (Figures 30 and 31) has an ascending growth curve and oscillating () due to the leading edge shock wave, which also causes an increase in the temperature, skin friction coefficient, and heat flux. The pressure continues to reduce until the beginning of the ramp. The ramp causes an increase in pressure and continues to increase until the end of the ramp. It is noticed that the increase in the Mach number and the angle of the ramp cause higher values of pressure on the wall.
The skin friction coefficient (Figures 32 and 33), the Stanton number (Figure 35), and the heat flux (Figure 34) reach the maximum value () when passing through the leading edge shock and continue decreasing along the plate; then, they have a slight increase on the ramp and finally maintain constant values until the end of the ramp.
The skin friction coefficient is lower in cases with an adiabatic wall (higher temperature). This is because the increase in the wall temperature increases viscosity which, in turn, increases hydrodynamic and thermal boundary layer thicknesses. Thus, the decrease in gradients close to the wall leads to a decrease in the skin friction coefficient. However, the thicker boundary layer is more susceptible to separation [14]. The increased Mach number freestream, compression corner angle, and adiabatic wall condition result in higher values of these properties.
The curves behave differently depending on boundary layer separation (Figures 36 and 37), which occurred in the cases of Mach 3 at 25 deg () and 30 deg () and Mach 4 (), 5 (), and 6 () at 30 deg for constant wall temperature conditions and Mach 3 (), 4 (), 5 (), and 6 () at 25 deg and Mach 3 (), 4 (), 5 (), and 6 () at 30 deg considering adiabatic wall conditions. That is, an increase in the Mach number, a decrease in the corner compression angle, and a decrease in the wall temperature interfere and reduce the possibility of separating the boundary layer. A more accentuated growth can be observed in the Stanton number and the skin friction coefficient in cases with boundary layer separation, see Figures 35–37, respectively.
Figures 38–45 show the pressure, temperature, Mach number, and numerical schlieren for Mach number 3 at 15 deg and the compression corner at constant wall temperature and adiabatic walls. It is possible to visualize the leading edge shock generated at the beginning of the ramp, with the deceleration of the flow and increase in temperature and pressure. There is also a subsonic laminar boundary layer (Figures 39, 40, and 44) on the ramp wall due to its no-slip condition.
In cases at 15 deg, the pressure gradient was not enough to cause boundary layer separation, and the ramp induces a pressure increase that causes compression waves to merge in a shock. Velocity (Figure 40) increases within the boundary layer until it reaches 99% freestream. This variation is greatest near the surface and varies more slowly at points further away. The shock becomes more intense as it moves away from the surface until reaching the inviscid solution. In numerical schlieren (Figures 41 and 45), it is possible to see the leading shock wave in both cases and the shock in the adiabatic wall (because the shock is stronger).
Comparing the pressure (Figures 38 and 42) and temperature (Figures 39 and 43) distributions, it can be seen that the cases with adiabatic walls have higher values than the cases with constant temperature. This happens because the adiabatic wall temperature is much higher than in the case where the wall temperature was taken as constant, . Thus, density has lower values, the boundary layer becomes thicker, the shock wave becomes stronger, and there is an increase in pressure within the shock layer [8].
Figures 46–53 show the pressure, temperature, Mach number, and numerical schlieren for the Mach number 6 at 15 deg and the compression corner at constant wall temperature and adiabatic wall. The physical behavior is the same for Mach 3 and 15 deg, and there is the formation of leading shock waves, compression waves, and shocks after the compression corner and the subsonic laminar boundary layer. However, in cases with Mach 6, there is a stronger leading edge shock (hypersonic), that is, because pressure and temperature values and compression waves are stronger. In cases with an adiabatic wall, it has higher pressure and temperature values than in cases with constant temperature.
Figures 54–61 show the pressure, temperature, Mach number, and numerical schlieren for Mach number 3 at 30 deg and the compression corner at constant wall temperature and adiabatic walls. It is possible to see the leading shock wave and the laminar boundary layer across the plate. In cases with 30 deg, the ramp causes a greater increase in the pressure gradient (Figures 54 and 58), causing the boundary layer to separate (S). In Figures 56 and 60, it is possible to observe the recirculation flow during the separation region and the reattachment point (R). The leading shock wave intersects with the reattachment shock in type VI shock-shock interference, causing a series of expansion waves to form.
This interference (Figures 54 and 58) takes longer to happen at higher Mach numbers due to the greater angle of the leading edge shock wave and in cases with constant wall temperature because the separation shock is weaker. In cases with constant wall temperature and adiabatic walls, one can also visualize the triple point T (Figures 54, 55, and 58). This point refers to the collision of the separation and reattachment shocks. It occurs at earlier stations in cases with constant wall temperature, and as it has a smaller separation region, reattachment also occurs earlier. With increasing wall temperature, the density (Figures 57 and 61) and Reynolds number decrease. As the Reynolds number decreases, the boundary layer resistance to the pressure gradient increases. Therefore, cases with constant wall temperature have a smaller separation bubble (Figures 54 and 58). The physical phenomena in the temperature and numerical schlieren distributions can be better visualized in the case of constant temperature.
Figures 62–69 show the pressure, temperature, Mach number, and numerical schlieren for Mach number 6 at 30 deg and the compression corner at constant wall temperature and adiabatic walls. The same phenomena are observed for Mach 3 and 30 deg: leading edge shocks, separation shocks, reattachment shocks, type VI shock-shock interference, laminar boundary layers, and expansion waves. In the case of Mach 6 with an adiabatic wall, one can also visualize the triple point T, that is, the interference of the separation shock with the leading shock wave and reattachment shock.
In Figures 64 and 68, it is possible to observe the recirculation flow during the separation region and the reattachment point (R). Increasing the Mach number and decreasing the temperature (in the case of constant temperature) reduce the temperature (Figures 63 and 67) and pressure gradient (Figures 62 and 66), make the shock region less intense, and decrease the separation region.
4. Conclusion
The objective of this paper was to study the limits of the geometry angles to predict the conditions in which the boundary layer separation occurs when it interacts with a compression corner. A code was developed using the MacCormack method to simulate the supersonic flow. The methodology and code for the transformation of variables were developed. Thus, simulations were performed for Mach numbers 3, 4, 5, and 6 and deflection angles of 10, 15, 20, 25, and 30 deg. Boundary layer separation was observed for Mach 3 at 25 deg and 30 deg and Mach 4, 5, and 6 at 30 deg for constant wall temperature conditions and Mach 3, 4, 5, and 6 at 25 deg and Mach 3, 4, 5, and 6 at 30 deg considering adiabatic wall conditions. It was observed that an increase in the Mach number, a decrease in the corner compression angle, and a decrease in the wall temperature interfere and reduce the possibility of separating the boundary layer. It was possible to visualize the leading shock wave and the laminar boundary layer across the plate. In cases with separation, the deflection of the ramp caused an increase in the pressure gradient, making it possible to visualize the boundary layer to separate, the recirculation flow region and the boundary layer reattachment. In some cases, the leading shock wave intersects with the reattachment shock in type VI shock-shock interference, causing a series of expansion waves to form.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this article.